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Abstract 


The  Vector  Slope  Gauge  (VSG)  is  a  3S-GHz  FM-CW  scatterometer  that  has  the 
unique  capability  of  simultaneously  (nearly)  measuring  the  range  and  backscattered 
power  to  three  points  on  the  ocean  surface.  With  three  ranges  and  knowledge  of  the 
experimental  geometry,  the  wave  height  at  each  footprint  can  be  obtained.  The  three 
footprints  form  a  plane  surface  which  enable  two  orthogonal  slope  components  to  be 
obtained.  Obtaining  the  slope  of  ocean  waves  is  important  because  it  is  correlated 
with  the  backscattered  power.  By  obtaining  the  vector  slope,  one  does  not  have  to 
make  any  assumptions  about  the  linearity  or  long-crestedness  of  the  ocean  waves. 

With  a  time  series  record  of  ocean  wave  heights  and  slopes,  one  can  learn  a  great 
deal  about  the  ocean  surface.  Spectral  analysis  of  the  recorded  time  series  yields 
information  about  the  wave  height  power  spectral  density,  the  mean  wave  direction 
vs.  fi«qtiency,  and  the  directional  width  spectrum.  T*  s  results  from  the  VSG  are 
similar  to  diose  obtained  fix>m  a  pitch*and-roll  buoy.  In  addition,  since  the  slope 
distribution  of  ocean  waves  is  nearly  normal,  the  moments  of  a  bivariate  normal 
distribution  can  be  used  to  fit  an  ellipse  to  the  wave  slopes.  The  orientation  of  the 
major  axis  of  the  ellipse  indicates  the  direction  of  dominant  wave  travel,  with  180** 
ambiguity.  The  ellipse  also  yields  information  about  the  statistics  of  the  slope 
distribution.  Due  to  the  asymmetry  of  the  waves,  the  center  of  the  ellipse  is  shifted 
fi»m  the  origin  in  the  direction  from  >^ch  the  dominant  waves  are  traveling.  Thus, 
the  wave  direction  ambiguity  can  resolved  by  fitting  an  ellipse  to  the  directional 
histogram  of  the  slope  distribution  using  the  least  square  method. 

As  with  any  instrument,  the  VSG  has  some  inherent  errors  due  to  the  method  of 
measurement  TlMse  errors  include  phase  shifts  in  the  recorded  time  series  due  to 
measuring  along  a  slant  range,  phase  shifts  due  to  non-simultaneous  measurements. 


and  errors  due  to  q)proxiinating  the  slope  at  a  point  with  a  plane.  Slant-range 
measurements  cause  the  measured  time  series  of  slopes  to  be  greater  than  or  equal  to 
the  actual  slope  at  all  times  in  any  direction.  The  overestimate  varies  from  0**  at  the 
mean  sea  level  to  2**  at  the  crest  and  trough.  The  derivative  approximation  has  little 
effect  on  either  the  wave  height  time  series  or  the  slope  time  series.  Non- 
simultaneous  measurements  cause  the  magnitude  of  the  slope  to  be  overestimated  in 
both  directions  at  all  times.  The  amount  of  overestimate  is  generally  around  2”.  This, 
although  small,  causes  an  error  in  the  determination  of  the  mean  wave  directiotL 
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1.  iBtrodaction 

1.1  Background  of  Saxon  experiment 

The  Synthetic  Aperture  radar  and  X>band  Ocean  Nonlinearities  (SAXON) 
e}q)erinient  took  place  in  November  of  1990  on  the  Nordsee  research  platform.  The 
tower  was  situated  in  about  30  meters  of  water  just  off  the  coast  of  Germany  in  the 
North  Sea.  The  University  of  Kansas  was  one  of  several  institutions  involved  in  the 
txjpensneat  The  purpose  of  the  experiment  was  to  study  microwave  backscatter  and 
SAR  images  of  the  ocean  sur&ce  for  high  sea  states  [Plant  &  Alpers,  1991].  To  that 
endeavor.  The  University  of  Kansas  operated  a  3S>GHz  scatterometer,  as  well  as  C- 
and  X-band  radars. 

The  3S-GHz  scatterometer  that  KU  operated  is  more  commonly  called  the 
Vector  Slope  Gauge  (VSG)  in  light  of  its  unique  ciq)abilities  (see  section  2.1).  The 


East 

Figure  1.1.  Diagram  of  Nordsee  platform  arui  location  of  VSG  and  10-m  anmnometer. 
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VSG  was  located  on  the  Northwest  comer  of  the  platform  about  20  meters  above  the 
mean  sea  level  (see  Figure  1.1).  The  VSG  was  pointed  in  the  general  direction  of 
oncoming  waves  at  moderate  angles  of  incidence.  To  avoid  cormption  of  the  data 
due  to  the  location  of  the  tower,  we  analyzed  data  collected  only  when  waves  were 
coming  between  225**  and  45°  (0°  is  North),  and  ^en  the  wind  was  within  the  same 
region.  The  wind  speed  was  obtained,  when  available,  fix>m  a  sonic  anemometer  S 
meters  above  the  mean  sea  surface  off  the  Northeast  comer  of  the  tower.  When  these 
data  were  unavailable  the  wind  speed  was  obtained  from  a  c\q)  anemometer  47  meters 
above  the  mean  sea  level  at  the  southwest  comer  of  the  platform.  An  algorithm  was 
then  used  to  estimate  the  wind  speed  at  5  meters. 

Ocean  measurements  were  also  made  with  a  Wavec  pitch-and-roll  buoy  operated 
by  die  Federal  Maritime  and  Hydrogrtqihic  Agency  of  Hamburg,  Germany.  The 
buoy,  located  near  the  tower,  recorded  time  series  of  heave  and  pitch  and  roll  angles. 

IJL  Importance  of  studying  ocean  waves 

In  recent  years,  concern  over  the  environment  has  increased  dramatically  with  the 
reveladrm  of  the  ozoae  "hole"  and  talk  of  tlw  greenhouse  effect  Since  70%  of  the 
eardi's  surfrue  is  covered  by  water,  understanding  the  air-sea  interfrue  plays  a  nujor 
role  in  developing  global  climate  mo<tels.  To  understand  how  wind  affects  ocean 
waves,  it  is  important  to  understand  how  die  energy  in  ocean  waves  is  distributed  with 
regard  to  frequency  and  direction  (the  directional  spectrum).  Knowledge  of  the 
directional  spectrum  of  ocean  waves  is  also  important  for  determining  the  effect  on 
flipping  [Neumann  &  Pierson,  1963],  on  off-shore  structures  [Kuik  et  al.,  1988], 
coa^  sedimmit  tran^rt,  and  wave  diffraction  and  refraction  from  the  shore 
[Trageser  &  Elwany,  1990]. 
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The  advantage  of  using  satellite  microwave  remote-sensing  instruments  to 
measure  ocean  parameters  is  that  they  can  be  used  in  any  type  of  weather.  The 
backscattered  power  fiom  the  radar  signal  is  affected  by  tilt  and  hydrodynamic 
modulation.  Thus,  to  better  understand  how  the  radar  signal  is  modulated  by  ocean 
waves.  The  University  of  Kansas  operated  several  radars  as  part  of  the  SAXON-FPN 
ejqperiment. 

13  Outline  of  thesis  and  summary  of  results 

The  objective  of  my  thesis  is  to  analyze  the  inherent  errors  involved  with  the 
VSG,  and  to  show  that  they  do  not  seriously  detract  from  the  VSG's  cq>ability  of 
determining  ocean  wave  heights  and  vector  slopes.  Errors  occur  with  the  VSG  due  to 
the  radar's  distance  measurement  along  a  slant  range,  approximation  of  the  derivative 
at  a  point  on  the  ocean  surfrax  by  the  slope  of  a  plane,  and  nonsimultaneity  of  the 
three  range  measurements.  These  errors  induce  phase  ^fts  in  the  measured  time 
series  of  wave  heights  and  slope.  Tlw  effect  on  the  slope  time  series  is  more 
significant  due  to  the  frwt  that  individual  range  errors  cannot  be  averaged. 

The  two-dimensional  distribution  of  wave  slopes  measured  by  the  VSG  is  nearly 
bivariate  Gaussian.  The  moments  of  the  slope  distribution  can  be  used  to  fit  an 
ellipse  to  the  wave  slopes.  The  mean  direction  of  wave  travel  is  given  by  tin 
(mentation  of  the  major  axis  of  the  ellipse,  with  180**  ambiguity.  Due  to  the 
asymmetry  of  ocean  waves,  however,  the  center  of  the  ellipse  is  shifted  finm  the 
(»rigin  in  tlm  direction  fiom  uhich  the  waves  are  coming.  Thus,  by  fitting  an  ellipse 
to  a  histogram  of  the  wave  slopes  the  directional  ambiguity  can  be  resolved. 

Reccuds  of  the  wave-height  time  series  and  slope  time  series  allow  calculation  of 
the  directional  spectrum  of  ocean  waves.  From  the  directional  spectrum,  the  wave 
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height  power  spectral  density,  the  mean  wave  direction  spectrum,  and  the  directional 
width  spectrum  can  all  be  calculated.  These  same  spectral  properties  are  routinely 
calculated  from  pitch-and-roU  buoy  data.  A  WAVEC  pitch-and-roll  buoy  was 
operated  by  the  Federal  Maritime  and  Hydrographic  Agency  of  Hamburg,  Germany 
during  the  SAXON-FPN  experiment.  The  data  from  this  buoy  were  analyzed  by  F. 
Ziemer  of  GKSS,  Geesthacht,  Germany,  and  the  results  were  made  available  to  us  for 
comparison.  The  spectral  parameters  determined  from  the  VSG  data  compare 
frivorably  with  those  from  the  pitch-and-roll  buoy. 

2.  Theoretical  Analysis  of  Errors 

2.1  The  Vector  Slope  Gauge 

The  VSG  is  a  35-GHz  FM-CW  scatterometer  that  uses  a  single  parabolic 
reflector  with  three  switchable  feeds  to  measure  simultaneously  the  range  to  and  the 


backscattered  power  from  three  closely  spaced  points  on  the  ocean  surface.  The  three 
feeds  are  situated  such  that  the  radar  footprints  form  a  right  angle  on  a  surface 
perpendicular  to  the  center  beam.  The  3-dB  beamwidth  of  each  beam  is 
iq^mximately  2**,  and  2.3**  separates  the  center  of  each  beam  from  the  vertex  of  the 
right  triangle.  The  VSG  was  mounted  roughly  20  meters  above  the  mean  sea  level. 
Thus,  at  47**  the  area  of  each  footprint  is  approximately  1.3  m^,  and  the  center-to- 
center  distance  between  footprints  in  each  leg  of  the  right  triangle  is  roughly  1 .8  m. 

The  beams  of  the  VSG  were  switched  at  30  Hz,  and  a  measurement  of  the 
intermediate  frequency  (IF)  and  backscattered  power  was  made  at  each  sampling 
interval.  Thus,  data  were  recorded  for  each  beam  at  10  Hz.  For  an  FM-CW  radar,  the 
difference  between  transmitted  frequency  and  received  frequency,  f),,  is  proportional 
to  the  range  to  the  target,  and  is  given  by 


wliere  B  is  the  bandwidth  of  the  sweep  signal,  r  is  the  range  to  the  target,  fg,  is  the 
modulation  rate,  and  c  is  the  speed  of  propagation.  The  VSG  was  swept  over  a  500 
MHz  bandwidth  centered  at  6  GHz.  The  modulation  rate  was  chosen  so  that  the  mean 
IF  frequency  was  455  kHz. 

With  the  range  to  three  points  on  the  ocean  surface,  and  knowledge  of  the 
orientation  of  the  radar,  the  height  of  the  ocean  surface  within  each  footprint  can  be 
calculated.  Since  the  three  points  define  a  plane,  the  slope  of  the  plane  can  be 
determined  with  regard  to  any  direction. 
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2J2  Coordinate  Transformations 

When  working  with  data  from  the  VSG,  one  needs  fiequently  to  transform 
between  different  coordinate  systems.  Figure  2.2.1  depicts  the  situation  when  the 
VSG  is  rotated  from  the  earth  coordinate  system  (non-primed  axes)  to  the  antenna 
coordinate  system  (double-primed  axes).  The  VSG  is  first  rotated  an  angle  of  9  about 


Figure  2.2.1  Transformation  of  radar  coordinate  system  to  earth  coordinate  system. 


the  X  axis,  and  then  an  angle  of  6  about  the  z*  axis.  To  determine  the  coordinate 
transformation  matrix  that  results  from  the  rotation  of  9,  see  Figure  2.2.2.  The 


Figure  2.2.2  Rotation  of  6  about  the  X  axis, 
coordinates  of  point  P  in  the  XYZ  coordinate  system  are: 


Px  =  0  /^  =  rcos(a)  /*r  =  rsin(a)  Eq  2.2.1 


In  the  X'YZ'  coordinate  system,  the  coordinates  of  point  P  are: 


Px'=Pk  /V'=rcos(a-9)  A’=rsin(a-0) 

/V'=  r(cos(a)cos(0) +sin(a)sin(0))  r(sin(a)cos(0)  -  sin(9)cos(a)) 
fycos(Q)  +  Asin(0)  Pk'=  Acos(0)  -  I)isin(Q) 


Thus, 


Eq.  2.2.2 
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To  determine  the  coordinate  transformation  matrix  that  results  fiom  the  rotation  of  5, 
see  Figure  2.2.3.  The  coordinates  of  point  P  in  the  X'VZ'  coordinate  system  are: 


Px'=r  cos(a)  Py=r  sin(a)  Pz'  =  0 

r 


Figure  2.2.3  Rotation  of  5  about  the  Z'  axis. 


In  the  X'Y'Z”  coordinate  system,  the  coordinates  of  point  P  are: 


’=rcos(a-8)  /y’=rsin(a-6) 

'=  r(cos(a)cos(5)  ■<-sin(a)sin(5))  Py-  r(sin(a)cos(5)  -sin(5)cos(a)) 

•*  PxcosiS)  /y  sin(5)  Py=  /y  cos(5)  -  PafsinCS) 

Thus, 


’A"' 

cos(5)  sin(5)  O' 
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A" 

= 

-sin(5)  cos(8)  0 

X 

Py' 

Pz" 

0  0  1 
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Eq.  2.2.4 


Pi"=  A’ 


Eq.  2.2.5 


Eq.  2.2.6 


Now,  to  transform  firom  the  earth  coordinate  system  to  the  radar  coordinate  system, 
multiply  the  coordinate  transformation  matrix  in  Eq.  2.2.6  by  the  coordinate 
tranrformation  matrix  in  Eq.  2.2.3. 
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Conversely,  if  the  coordinates  are  known  with  respect  to  the  radar,  then  the  matrix  to 
transform  them  to  the  earth  coordinate  system  is: 


Px' 

cos(8) 

-sin(8)  0 

Px" 

Py 

sin(5)cos(0) 

cos(8)cos(0)  -sin(0) 

X 

Py" 

Eq.  2.2.8 
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23  Measurement  of  Wave  Height  and  Slope 

The  vector  slope  gauge  (VSG)  measures  the  range  to  three  separate  points  on  the 
ocean  surface  through  the  use  of  its  switched-beam  antenna.  Given  the  ranges,  the 
location  of  each  point  in  the  antenna  coordinate  system  is  given  by  (see  Figure  2.2.1): 

A  =  (/^sinOi2)cos(90-y),  r^sin(Pi2)sin(90-y),  -r,cos(Pi2)  ) 

B  =  (0,  0,  -z^)  Eq.  2.3.1 

C  *  (  0,  Z5  sin(P23  ^  cos(P23  )  ) 
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To  determine  the  wave  heights  with  respect  to  the  mean  sea  level,  the  coordinates 
in  Eq.  2.3.1  are  first  transformed  to  the  earth  coordinate  system  using  Eq.  2.2.8. 
Then,  since  the  radar  is  located  approximately  20  meters  above  the  mean  sea  level,  20 
meters  is  added  to  the  z  (earth)  coordinate  of  each  point. 

Since  the  locations  of  three  points  on  the  ocean  sur&ce  are  known,  two 
orthogonal  slope  components  can  be  obtained.  The  normal  to  the  plane  containing 
the  three  points  is  given  by  (see  Figure  2.2.1) 


N  =  BA  y^BC 


Eq.  2.3.2 


The  slope  of  the  plane  in  the  y-direction  is  tan"^  ( — —)  and  the  slope  of  the  plane  in 


N, 


I  -Nr 

the  x-direction  is  tan  ( — ^). 

N, 


The  Matlab  program  slopeab.m  was  written  by  Sam  Haimov  and  Vahid  Hesatry 
to  determine  orthogonal  components  of  the  slope,  given  the  range  measurements.  In 
slopeab.m,  the  coordiruUe  transformation  used  is  different  than  the  one  given  in  Eq. 
2.2.7  or  2.2.8.  The  two  coordinate  transformation  matrices  give  the  same  results, 
however,  because  the  angles  6  and  9  in  slopeab.m  (corresponding  to  9  and  5  in  the 
derivation  above)  ate  positive  when  measured  clockwise  in  program  slopeab.m  as 
opposed  to  counter-clockwise  in  the  above  derivation. 


2.4  Inherent  errors  in  VSG  measurements 

In  every  instrument  there  are  measuremrat  errors  due  to  noise,  and  potoitial 
mors  due  to  calibration  problems.  With  the  VSG,  there  ate  also  errors  that  occur 
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completely  independent  of  noise  and  calibration  problems.  These  errors  are  due 
simply  to  the  way  the  measurements  are  made  and  some  of  them  cannot,  in  general, 
be  reduced  purely  by  system  design.  (It  may  be  possible,  however,  to  reduce  them  by 
changing  the  experimental  geometry.)  The  inherent  errors  are  due  to  measuring  along 
a  slant  range,  approximating  the  slope  at  a  point  with  a  plane,  and  non-simultaneous 
measurements. 

In  the  ideal  situation  and  neglecting  noise,  the  measured  time  series  of  a  pure 
sinusoidal  wave  of  one  frequency,  direction,  and  amplitude  mcident  iq)on  the  radar 
measurement  site  would  be  an  exact  replica  of  the  actual  wave.  Due  to  tl^  inherent 
errors,  however,  phase  shifts  occur  in  the  measured  time  series  of  wave  heights  (see 
Figure  2.4.1).  The  measured  wave  l^ight  time  series  depicted  in  Fig.  2.4.1  was 
calculated  from  the  MATLAB  program  range.m  and  can  be  found  in  Appendix  A.  It 

Maasurad  vs.  Actual  Wave  Haight  Tima  Saiias 
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Figure  2.4.1.  Exanq>le  of  phase  error  in  3*beam  average  wave  height  time  series  for  a 
steep  wave  due  to  slant-range,  non-simultaneous  measurement.  Incidence  angle,  47**; 
Antoina  hei^t,  20  m;  Wave  frequency,  0.2  Hz,  6  -  0**,  Upwave  look  direction. 
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Msasured  vs.  Actnl  Slops  Tims  Ssriss 
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Figure  2.4.2.  Example  of  total  error  in  slope  time  series  due  to  slant-range,  ixm- 
simultaneous  measurement,  and  the  derivative  {q)proximation.  Incidence  angle,  47**; 
Anteruui  height,  20  m;  Wave  frequency,  0.2  Hz;  Wave  height,  3  m;  5  »  0**;  Upwave 
look  direction. 
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takes  into  account  measurement  along  the  slant  range  as  well  as  non-simultaneous 
measurements.  The  maximum  phase  occurs  at  the  crest  and  the  trough  because 
at  these  points  the  measurement  site  is  a  maximum  horizontal  distance  fiom  the  tnw 
location. 

In  tlM  slope  time  series,  the  individual  errors  in  each  range  measurement  have  a 
greater  effect  because  they  are  not  averaged.  Measuring  along  the  slant  range  induces 
a  i^iase  shift  in  the  recorded  time  series;  approximating  the  derivative  causes  die 
recorded  time  series  to  be  modulated  by  a  sine  function;  and  ncm-simultaneous 
measurements  cause  the  estimated  wave  direction  to  be  in  error.  The  total  effect  of 
these  oTors  is  depicted  in  Figure  2.4.2.  Since  the  oncoming  wave  is  headed  directly 
toward  the  radar,  diete  is  no  slope  in  the  cross  (x)  direction.  Note  that  for  this 
particular  exanq>le,  however,  the  measured  slope  in  the  x  direction  can  be  as  much  as 
±3^  As  a  result,  the  mean  wave  direction  would  be  recorded  as  coming  fiom 

90*-tan"'^yj»ll® 

an  error  of  1 1^.  Now  that  the  total  error  can  be  seen,  it  is  worth  looking  at  eadi 
source  of  error. 

2,5  Error  in  wave  hc^t  time  scries  due  to  slant-range  measnrement 

A  radar  measures  the  distance  to  a  particular  point  along  the  radar  beartL  Thus, 
the  VSG  measures  the  distance  to  the  pointi  of  intersection  between  eadi  of  the  radar 
beams  and  the  ocean  surfice.  Since  the  ocean  surfice  moves  due  to  the  waves,  the 
points  of  measurement  move  along  the  radar  beam  rather  Aan  vertically.  This  causes 
a  phase  diift  in  the  time  series  of  ranges  which  affects  the  determination  of  ocean 


wave  heights  and  slope.  To  estimate  the  amount  of  enor  in  the  wave  height  time 
series  and  slope  time  series  induced  by  the  non-vertical  displacement  for  waves  of 
various  frequencies  and  directions,  it  is  necessary  to  simulate  the  ocean  surface 
Simplifying  a  great  deal,  the  surfrkce  of  the  ocean  can  be  (tescribed  by: 

z  =  i4cos(2i^'i-i:(sin^)x  +  ib(cos^)y)  Eq.  2.5.1 

where  z  is  the  wave  height,  A  is  the  wave  amplitude,  f  is  the  wave  frequency,  k  is  the 
wave  number,  and  ^  is  the  angle  (measured  clockwise)  between  the  +y-axis  (RLD) 
and  the  directimi  from  which  the  waves  are  coming.  If  we  assume  that  the  origin  of 
the  earth  coordinate  system  is  located  where  the  radar  signal  originates  with  +z 
indicating  an  up  direction,  then  the  z  coordinate  of  any  point  on  the  ocean  surfrtte  is 

z*i4cos(2i?/f+k(sin^)*  +  k(cos4)y)-A  Eq.  2.5.2 

udiere  h  is  the  height  of  the  antenna  above  die  mean  ocean  surfiwe.  Now,  fixim  Eq. 
2.2.8  and  Eq.  2.3.1,  die  locations  of  die  points  of  measurement  in  die  earth  coordinate 
system  are  given  by  Eq.  2.5.3.  The  ai^noimate  x,  y,  and  z  values  from  Eq.  2.5.3  can 

*1  *  n  [«>a(S)sin(Pi2  )co9(r  -  90) + sin(5)stnO|2  )siii(y  -  90)] 
yi  «n[coa(e)sin(5)siii(p,2)cos(y  -90)-co8(5)cos(e)sin(Pi2)sin(T  -90)+sin(e)coa(P,2)] 
zj  =ii[siii(e)siii(5)8inOi2)cos(y  -9O)-cos(5)siii(0)suiO,2)sin(y  -90)-cos(e)cosO,2)] 

X2  =0 

since) 

X3  »/5(sin(5)sin023)] 

3^3  *r5[cos(8)cos(O)sin(p23)  +  sin(0)cos(p23)l 

Z3  « [cos(5) since)  sinCP23 )  ^  cos(0)cos(P23  )]  2.5.3 
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be  substituted  into  Eq.  2.5.2,  resulting  in  three  equations  and  three  unknowns. 
Newton's  mediod  is  then  used  to  solve  for  the  three  unknown  raises.  Thus, 


'ii+l 


Eq.  2.5.4 


where  is  the  n^  value  of  the  range. 

In  the  sinq>le  case  of  waves  of  one  fiequency  at  moderate  angles  of  incidence, 
(me  does  not  have  to  worry  about  local  minima  vs.  global  minima  because  each  radar 
beam  intercepts  the  ocean  surface  at  only  one  point.  However,  if  waves  of  multiple 
frequencies  are  added  together,  or  if  the  angle  of  mcidence  is  close  to  grazing  angles, 
one  needs  to  check  to  make  sure  diat  dw  range  given  as  die  soludim  to  the  problem  is 
the  actual  range  from  Mhich  the  radar  signal  came.  The  correct  range  should  be  the 
shcntest  range  i^ch  is  a  soluticm  to  the  nonlinear  equadcm  given  above. 

The  program  slrange.m  was  written  in  MATLAB  3.5  to  determine  die  x,  y,  and  z 
eardi  coordinates  of  die  points  of  intersection  and  also  die  ranges  to  the  points  of 
intersectitm.  It  assumes  diat  each  set  of  three  measurements  is  made  simultaneously. 
To  determine  die  fdiase  error  in  wave  Iwight  measurements  associated  widi  the  non- 
vertical  diqilacement  of  the  point  of  measurement,  one  needs  to  determine  the  wave 
heights,  assuming  diat  die  point  of  iiKasurement  is  moving  vertically.  The  MATLAB 
3.5  pst^ram  vrange.m  was  written  to  determine  the  x,  y,  and  z  earth  coordinates  of 
the  points  of  intersection  and  also  the  ranges  to  the  points  of  intersection  givm  the  x 
and  y  coordinates  olrtained  from  slrange.m  and  several  odier  variables.  The  x  and  y 
coordiiuttes  are  required  so  diat  the  data  from  die  two  programs  can  be  properly  (diase 
matched  (see  Figure  2.5.1).  That  is,  the  range  time  series  and  the  wave  height  time 
series  of  die  dnee  points  will  be  in  phase  exacdy  at  the  mean  sea  level. 


Figure  2.S.1  Mean  time  series  of  simulated  wave  height  measuranents  illustrating 
the  jriiase  error  due  cmly  to  a  non-vertically  moving  point  of  measurement  The  solid 
line  represents  the  average  of  each  measuraimnt  made  along  the  radar  beam.  The 
dashed  line  represents  die  average  of  three  vertical  beam  measurements.  Incidence 
angle,  47*;  Antenna  height  20  m;  Wave  fiequency,  0.1953125  Hz,  5  «  0*,  Upwave 
look  direction. 


As  previously  noted,  phase  shifts  in  die  time  domain  result  in  harmonics 
appearii^  in  the  fiequency  domain  (see  Figure  2.5.2).  The  particular  example  shown 
in  Fig.  2.5.2  is  for  the  same  wave  height  time  series  depicted  in  Fig.  2.5.1.  It  is  for  a 
particulariy  steep  wave;  the  wave  height  to  wave  length  ratio  in  this  case  is  near  die 
theoretical  limit  of  1/7.  As  can  be  seen,  even  for  such  an  extreme  case  die  second 
harmonic  is  small,  and  higher  order  harmonics  are  even  smaller. 
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Figure  2.S.2.  Harmonics  occuning  as  a  result  of  slant-range  measurement  Incidence 
angle,  47**;  Antenna  height,  20  m;  Wave  frequency,  0.19S312S  Hz;  Wave  height,  3  m; 
5*0*,  Upwave  look  direction. 
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The  maximum  phase  error  for  a  single  beam  >^ch  occurs  at  the  crest  and  trough 
can  be  determined  from  the  amount  of  horizontal  shift  which  occurs  in  the  radar  look 
direction  as  a  result  of  slant-range  measurement.  According  to  Figure  2.S.3,  the 
horizontal  phase  shift  along  the  direction  of  wave  propagation  is  i4tan(6).  The 


Figure  2.5.3.  Geometry  used  to  calculate  maximum  amount  of  phase  error  in  a  single 
beam  occurring  at  the  crest  and  trough  due  to  slant-range  measurement. 


component  of  diat  phase  shift  lying  in  the  radar  look  direction  is  simply 
i4tan(9)cos(^).  Thus,  the  maximum  amount  of  phase  error  for  a  single  beam,  in 
degrees,  is 


i4tan(9)cos(4)36y 

X 


Eq.  2.5.5 


When  the  wave  heights  from  the  three  beams  are  averaged  together  to  form  the 
mean  wave  height  time  series,  the  maximum  phase  error  for  the  mean  wave  height 
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time  series  will  be  slightly  different  than  this.  It  will  generally  be  less  than  the  error 
in  a  single  beam  unless  the  wave  direction  is  perpendicular  to  the  radar  look  direction. 

The  examples  shown  in  Figures  2.5.1  and  2.S.2  are  worst  case  scenarios  for 
several  reasons.  The  radar  is  looking  directly  at  the  oncoming  waves  so  that  the 
cos(^  +  term  is  1 .  Both  figures  are  also  for  steep  waves;  for  the  givm  X,  A  is  near 
its  maximum  theoretical  limit.  From  Eq.  2.2.5,  one  can  see  that  for  a  wave  of  a  given 
length  and  amplitude  and  for  a  given  incidence  angle,  the  phase  error  will  be  greater 
for  a  wave  in  the  vq>wave  or  downwave  look  direction  dian  for  any  other  direction. 
For  waves  coming  from  other  directions,  the  component  of  horizontal  shift  in  the 
radar  look  direction  is  less  than  the  horizontal  shift  in  the  direction  of  wave 
propagation.  For  a  wave  of  a  given  frequency,  amplitude,  and  direction,  the  phase 
error  will  be  worse  for  a  larger  angle  of  incidence.  As  the  incidence  angle  increases, 
Aere  is  more  horizontal  shift  in  the  radar  look  direction. 


Fraquancy  (Hz) 


Afnpiilude  (m) 
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Figure  2.5.4.  Maximum  phase  error  in  wave  iMight  time  series  due  to  slant-range 
measurement  Incidence  angle,  47**;  Antenna  height  20  m;  5  -  0**;  Upwave  lo(^ 
directioiL 
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Table  2.5.1  lists  the  maximum  phase  error  for  waves  of  various  amplitudes  and 
frequencies  given  that  the  angle  of  incidence  is  47**  and  the  direction  of  wave 


propagation  is  frx>m  0**.  The  values  are  also  plotted  in  Figure  2.5.4. 

Table  2.5.1.  Maximum  phase  error  (deg)  due  to  slant-range  measurements  for  waves 
of  various  frequencies  and  amplitudes.  Incidence  angle,  47**;  Anteima  height,  20  m; 
8  =  0*;  Upwave  look  direction. 


Amplitude  of  Ocean  Waves  (meters) 


X(m) 


243.7 


156.0 


108.3 


79.6 


f(Hz) 


Amplitude  of  Ocean  Waves  (meters) 


X(m) 


243.7 


156.0 


108.3 


79.6 


f(Hz) 


0.08 


.10 


0.12 


0.14 


0.16 


0.18 


0.20 


0.22 


2.6  Error  in  wave  height  time  series  due  to  non-simultaneous  measurement 


The  VSG  records  the  range  and  backscattered  power  from  the  ocean  surface 
every  33  ms.  Thus,  during  the  time  between  measurements  the  surface  moves 
slightly.  This  movement  shows  up  as  an  additional  phase  shift  in  the  measurements 
from  beams  2  and  3.  For  example,  if  the  waves  are  coming  directly  at  the  radar, 
height  measurements  from  beam  I  and  2  should  be  identical.  Due  to  the  finite 
switching  time  between  feeds,  however,  the  wave  height  measured  by  beam  2  will  be 
greater  than  the  wave  height  measured  by  beam  1  along  the  front  face  of  the  wave, 
while  along  the  back  face  of  the  wave,  the  wave  height  measured  by  beam  2  will  be 
less  than  that  measiuod  by  beam  1. 

By  comparing  Fig.  2.5.1  with  Fig.  2.4.1,  one  can  see  that  most  of  the  error  in  the 
wave  height  time  series  is  a  result  of  measuring  along  a  slant  range.  The  error  in  the 


Figure  2.6.1.  Difference  between  mean  wave  height  time  series  measured 
simultaneously  along  a  slant  range  and  measured  non>simultaneously  along  a  slant 
range.  The  sampling  interval  of  the  non-simultaneous  measurements  is  33  ms. 
IiKid«Dce  angle,  47**;  Antmna  height,  20  m;  Wave  frequency,  0.2  Hz;  Wave  height  = 
3  m;  8  *  0®;  Upwave  look  direction. 
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mean  wave  height  time  series  resulting  from  non-simultaneous  measurement  is  not 
very  significant  because  individual  beam  errors  ate  averaged  out.  Figure  2.6.1  shou^ 
the  difference  in  the  mean  wave  height  time  series  between  measurements  made 
simultaneously  along  a  slant  range  and  measurements  made  non-simultaneously  along 
a  slant  range.  Once  again,  this  is  for  a  steep  wave.  For  a  wave  having  less  curvature, 
the  difference  will  be  even  less.  Since  the  errors  in  the  mean  wave  height  time  series 
due  to  non-simultaneous  measurements  are  so  small,  there  is  no  point  in  investigating 
it  any  further. 


2.7  Error  in  Slope  Time  Series  due  to  Slant  Range  Measurement 
Measuring  along  a  slant  range  causes  errors  in  the  slope  time  series  similar  to  those 
that  occurred  in  the  wave-height  time  series.  Once  again,  there  is  a  phase  shift  in  the 
recorded  time  series  (see  Figure  2.7.1).  In  the  slope  time  series,  however,  the 
maximum  phase  shift  occurs  \riien  the  measured  slope  is  0**.  The  slope  is  0  nlien  the 
measuring  site  is  at  the  crest  or  trough  of  the  wave;  the  point  where  maximum  phase 
error  occurred  in  the  wave  height  time  smes.  When  the  measuring  site  is  at  the  mean 
sea  level,  the  phase  error  in  the  slope  time  soies  is  0  just  as  it  is  in  the  wave  height 
time  series. 

Another  important  observation  from  Fig.  2.7.1  is  fiiat  the  measured  slope  is 
always  greater  than  or  equal  to  the  actual  slope.  The  amount  of  difference  can  be  seen 
in  Fig.  2.7.2.  During  the  time  when  the  wave  is  above  the  mean  sea  level,  the  actual 
measurement  site  is  closer  to  the  radar  than  it  would  be  if  one  were  measuring 
vertically.  The  slope  of  the  positive  portion  of  a  sinusoidal  wave  is  a  continuously 
decreasing  function.  Thus,  since  the  actual  measurement  site  is  closer  to  the  radar. 
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Figure  2.7.1.  Error  in  slope  time  series  diw  only  to  slant-range  measurement  Both 
time  series  were  obtained  fiom  simultaneous  measurements,  and  slopes  were 
calculated  from  die  plane  formed  by  the  foo^nrints.  Incidence  angle  for  slant  range 
measurement  47*;  Antenna  height  20  m;  Wave  frequency,  0.2  Hz;  Wave  height  =  3 
m;  5  »  0^;  Upwave  look  direction. 

the  measured  slope  is  greater  than  it  should  be.  When  the  wave  is  below  tte  mean  sea 
level,  the  measurement  site  is  frffther  fiom  the  radar  than  it  would  be  if  one  were 
measuring  vertically.  The  slope  of  the  negative  portion  of  a  sinusoidal  wave  is  a 
continuously  increasing  function.  Thus,  since  the  actual  measurement  site  is  ferther 
fiom  the  radar,  the  measured  slope  is  once  again  greater  than  it  should  be.  It  is  also 
interesting  to  note  fiom  Fig.  2.7.2  that  the  difference  in  slope  between  slant  range 
measurements  and  vertical  measurements  occurs  at  the  2°4  harmonic  of  the  fipequency 
of  the  wave. 

Equation  2.S.S  gives  the  amount  of  phase  error  in  one  height  measurement  for  a 
sinusoidal  wave  of  any  length,  amplitude,  and  direction.  Obtaining  a  similar 
e)qnession  for  the  mean  wave  height  time  series  or  slope  time  series  is  not  practical 
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Figure  2.7.2.  Difference  in  slope  between  slant  range  measurement  and  vertical 
measurement  for  the  time  series  depicted  in  Fig.  2.7.1. 

because  of  the  interaction  of  the  three  beams  with  the  ocean  surface.  However,  Eq. 
2.S.S  can  be  used  qualitatively  to  estimate  how  much  phase  error  will  occur  in  the 
slope  time  series.  One  can  see  from  the  equation  that  larger  incidence  angles  and 
steeper  waves  (higher  A  to  X,  ratios)  will  result  in  more  phase  error. 

Both  slope  time  series  shown  in  Fig.  2.7.1  were  calculated  from  the  plane  formed 
by  the  duee  radar  footprints.  The  MATLAB  3.5  program,  slopeab.m,  determines  the 
slope  time  series  given  the  three  time  series  of  ranges,  the  incidence  angle,  and  the 
angle  of  rotation  about  the  z'  axis.  However,  \^en  all  three  beams  are  assumed  to 
measure  vertically,  as  in  Fig.  2.7.1  and  Fig.  2.7.2,  the  x  and  y  coordinates  of  each 
point  of  measurement  are  fixed.  (By  setting  the  incidence  angle  to  0^  only  beam  2  is 
truly  vertical.  The  x  and  y  coordinates  of  beams  1  and  3  will  change,  albeit  very 
slightly.)  Thus,  the  incidence  angle  would  have  to  continuously  change  in  order  for 
die  geometry  to  be  correct.  To  avoid  this  problem,  the  program  slopeab2m  was 
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Figure  2.7.3.  Error  in  slope  time  series  due  only  to  slant-range  measurement  Both 
time  series  were  obtained  fiom  simultaneous  measurements.  Incidence  angle  for  slant 
range  measurement,  47**;  Antenna  height  20  m;  Wave  frequency,  0.2  Hz,  Wave 
height  -  3  m,  6  *  0*,  Wave  direction,  30«,  RLD,  0*. 
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Mrritten  and  is  given  in  Appendix  A.  Slopeab2.m  calculates  the  slope  time  series 
givoi  the  X,  y  and  z  coordinates  of  the  points  of  intersection. 

In  Fig.  2.7.1,  only  the  slope  in  the  y  direction  is  shown.  In  that  example,  the 
waves  were  coming  directly  at  the  radar  so  there  was  no  slope  in  the  ooss  direction. 
In  Fig.  2.7.3,  the  waves  are  coming  from  30**,  and  all  other  characteristics  of  the 
waves  are  the  same  as  in  Fig.  2.7.1.  As  a  result,  there  is  slope  in  the  x  direction.  As 
can  be  seen  in  the  figure,  the  phase  error  in  the  x  direction  is  similar  to  the  error  in  die 
y  direction.  Thus,  although  the  VSG  accurately  measures  the  slope  (neglecting  errors 
due  to  die  derivative  r^iinoximation  and  non-simultaneous  measurements),  the  phase 
shift  due  to  slant  rai^e  measurements  causes  the  measured  slope  to  be  greater  than  the 
actual  slope  at  every  point  except  at  the  mean  sea  level. 

2^  Error  in  Slope  Time  Scries  due  to  First  Derivative  Approximation 

&TOT  occurs  in  the  slope  time  series  i^wnever  die  plane  used  to  calculate  the 
slope  is  tangent  to  a  point  other  than  the  one  for  i^ch  the  slope  is  being  estimated. 
The  estimate  of  the  slope  in  any  given  direction  is  [Hesany,  1994] 

Ar  Ar 

z(r+— ,r)-z(r-— ,0 

s(r,0  = - 2 - ^ 1 -  Eq.  2.8.1 

or 

assuming  the  plane  is  tangent  to  the  mic^int  of  the  radar  footprints.  The  true  wave 
height  is  given  by 

z(r,t)  =  Acos(jcat  +  kr)  Eq.  2.8.2 
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from  ^ch  tbe  slope  is  obtained  as: 


s(r,t)  =  -k4sin(ait  +  kr)  Eq.  2.8.3 

Substituting  Eq.  2.8.2  into  Eq.  2.8.1  results  in: 


Figure  2.8.1.  Location  of  point  P  vdiich  gives  the  minimum  error  in  the  derivative 
iqiproximation. 
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Equation  2.8.1  assumes  that  we  are  estimating  the  slope  at  the  midpoint  of  the 
three  radar  fooqnints.  Actually,  the  point  within  the  duee  radar  fooqnrints  tangent  to 
the  plane  is  not  stationary  with  respect  to  the  three  footprints.  Of  course,  if  it  were 
stati<»ary,  there  would  be  no  error  in  the  first  order  q>pfOximation  of  the  derivative. 

The  location  of  the  point  vibett  the  slope  is  being  estimated  that  gives  the 
minimum  amount  of  error  in  the  first  ordo*  aj^voximation  of  the  derivative  depends 
iqxm  tte  direction  of  oncoming  waves  with  respect  to  the  radar  look  directiotL  The 
q)|Mroximation  will  have  a  minimum  amount  of  oror  for  slopes  with  req)ect  to  the 
RLD  if  die  point  lies  along  a  line  perpendicular  to  tte  wave  direction  and  passing 
duou^  point  D  (see  Figure  2.8.1).  Conversely,  the  ap^Hoximation  will  have  a 
minimum  amount  of  ernxr  for  slopes  widi  reflect  to  the  cross  directitm  if  the  pomt 
lies  along  a  line  perpendicular  to  the  wave  direction  and  passing  duou^  point  E.  The 
point  of  minimal  error  for  both  time  series,  then,  lies  alrmg  a  line  midway  between 
and  parallel  to  the  lines  passing  through  D  and  E.  All  such  lines  will  pass  dirough 
point  P.  Howevo*,  if  the  radar  is  looking  directly  at  die  oncoming  waves  then,  for  die 
unidirecticmal  waves  bdng  considered  thus  &r,  there  is  no  slope  in  the  cross  directioiL 
Fot  such  a  case,  point  D  will  give  the  minimum  error  in  the  first  order  iqiproximation. 

If  measurements  are  made  with  all  diree  beams  vertical,  Ar  in  Eq.  2.8.4  is  a 
constant  In  such  a  case  the  error  in  die  first  order  rqiproximation  dqiends  only  on  the 
fi«quency  of  the  waves  and  Ar.  During  the  SAXON-FPN  experiment,  a  cormnon 
distance  between  die  radar  footprints  at  the  mean  sea  level  was  1.8  m.  Figure  2.8.2 
shows  the  relative  error  (  f(r,r)- j(r,/)]/5(r,r)  )  in  the  first  order  iqqiroximation 
of  the  derivative  for  waves  of  various  fipequendes. 


Derivativtt  Error  for  Vortical  Maaouramants 


Figure  2.t2.  Relative  errcv  in  the  first  order  approximation  of  die  derivative 
assuming  all  three  beams  are  vertical.  Ar»  1.8  meters. 

AldKM^  Fig  2.8.2  shows  the  error  awuming  all  three  beams  ate  vertical,  actual 
measuranents  were  made  along  a  slant  range.  Thus,  Ar  contmuomly  changes  and  is  a 
fimctkm  of  die  fduse  of  die  ocean  wave. 

Figure  2.8.3  shows  the  derivative  error  fiv  waves  of  various  frequencies  and 
heights  assuming  die  measuranaits  were  made  almig  a  slant  range.  The  higher 
frequency  waves  are  steqi-dieir  cotre^nding  heights  ate  large  relative  to  their 
Imgdi.  The  lower  frequmcy  waves,  vdiile  not  steep,  are  near  die  maximum  height 
opected  of  waves  near  die  Ncndsee  research  platform. 

As  previously  staled,  Ar  varies  and  u  a  fimction  of  the  jdiase  of  die  ocean  wave. 
Figures  2.8.4a,  b,  and  c  illustrate  die  location  and  the  amount  of  the  errcv  relative  to 
die  position  almig  the  wave.  Figure  2.8.4a  is  die  wave  height  time  series  frv  a  steqi 
0.4HzwaveMdientheradarisloddnginthetq>wavedirecti<HL  The  time  soies  from 
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Dcrivativ*  E  rror  W«v«  Ff.q.  (Hz). 


Wav*  Haight  (m) 


Figure  2.8.3.  Eiror  in  the  first  of  order  q)proximation  of  the  derivative. 
Measurements  woe  made  along  a  slant  range.  Incidence  angle,  47";  Antenna  height, 
20  m;  8  *  0*,  Wave  dire^mi,  0®,  RLD,  0®. 


beam  2  is  not  shown  because  it  is  identical  to  die  time  series  fcom  beam  1.  Figure 
2.8.4b  shows  the  corre^Kmding  slope  time  series  from  the  plane  and  fiom  the  point  D 
in  Fig.  2.8.1. 

From  Fig.  2.8.4,  mie  can  see  that  the  largest  error  in  the  derivi^ve  qiproximation 
occurs  near  die  trough  of  the  wave.  At  this  location  Ar  is  the  greatest  because  the 
fix^Mfints  are  a  maximum  distance  fiom  the  radar.  The  smallest  error  occurs  near  the 
crest  of  the  wave  because  at  this  location  the  foo^mnts  are  a  minimum  distance  from 
the  radar. 

From  these  results,  one  can  see  that  the  error  due  to  die  first  order  iqiproximation 
is  not  very  significant  Even  for  the  steepest  wave,  the  maximum  error  of  9%  is  less 
dian2®.  Furdmmore,  the  majority  of  energy-bearing  waves  during  theSAXON-FPN 
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0  36  72  108  144  180  216  252  268  324 


Phasa  (dag} 


Figure  2.8.4.  Wave  height  time  series  (a),  corresponding  slope  time  smes  (b),  and 
relative  error  due  to  derivative  qjproximation  (c).  Measurments  were  made  along  a 
slant  range.  Incidence  angle,  47*;  Antenna  height,  20  m;  5  »  0*,  Wave  frequency,  0.4 
Hz,  Wave  direction,  0*,  RLD,  0*. 


Figure  2.8.4b.  Slope  time  series  corresponding  to  the  wave  height  time  series  in  (a). 
This  iUustrates  udiere  along  the  wave  the  error  is  located. 
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Figure  2.8.4c.  Relative  error  in  derivative  q>proxiination  corresponding  to  the  wave 
height  time  series  in  (a).  This  illustrates  vdim  along  the  wave  tte  error  is  located. 


e}q)erunent  were  between  0.1  Hz  and  0.2  Hz.  For  waves  in  this  frequency  band,  the 
derivative  error  is  less  than  1%. 


2.9  Error  in  Slope  Time  Series  doe  to  Non-Simnltaneous  Measurements 

Errors  also  occur  in  the  slope  measurements  due  to  the  fact  that  the  range 
measurements  are  not  simultaneous  but  instead  are  sampled  33  ms  iq)art.  Although 
33  ms  is  a  small  amount  of  time  compared  to  the  period  of  waves  which  are  measured 
(2.5  sec  -  10  sec),  enough  movement  in  the  ocean  surface  occurs  during  this  time  so 
that  this  problem  cannot  be  neglected.  The  finite  switching  time  between  feeds 
causes  a  phase  drift  in  the  recorded  time  series  of  beams  2  and  3  (recall  section  2.6). 
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Since  the  delay  between  the  measurement  at  beam  3  and  the  measurement  at  beam  1 
is  66  ms,  the  phase  error  of  beam  3  will  be  twice  as  great  as  the  oror  in  beam  1  (see 
Figure  2.9.1). 


Phase  Difference  due  to  Finite  Sampling  Time 


Frequency  (Hz) 


Figure  2.9.1.  Phase  difference  in  ocean  wave  due  to  a  finite  sampling  time.  This  data 
assumes  that  the  footprints  are  fixed  with  respect  to  one  another. 

The  error  due  to  non-simultaneous  measurements  depends  iqwn  the  wave 
steepness.  For  long,  low  fiequency  waves,  the  height  does  not  change  as  rapidly  as  it 
does  for  short,  high  fiequency  waves.  However,  the  real  importance  here  isn't  the 
relative  error  but  the  absolute  error.  In  section  2.6,  it  was  noted  that  the  finite 
sampling  time  did  not  have  much  of  an  effi^  on  the  mean  wave  height  time  series. 
From  Figure  2.9.2,  one  can  see  that  for  these  steep  waves  the  error  due  to  a  finite 
sampling  time  is  on  the  order  of  several  centimeters.  (Note  that  the  errors  shown  ate 
only  for  a  sampling  interval  of  33  ms.  Beam  3  will  have  twice  as  much  error  shown 
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because  it  is  sampled  66  ms  after  beam  1.)  This  is  not  very  significant  for  the 
calculation  of  the  mean  'M^ve  height  time  series,  but  it  is  significant  for  the  calculation 
of  the  slope  time  series.  The  absolute  error  is  greater  for  the  low  fivquency  waves 
because  they  have  higher  wave  heights. 


Wave  Height  Error  due  to  Finite  Sampiing  Time 
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Figure  2.9.2.  Absolute  error  in  measured  wave  heights  for  waves  of  various 
fiequencies  and  heights.  The  error  is  due  to  a  finite  sampling  interval.  Beams  1  and  2 
are  assumed  to  measure  vertically. 


If  the  footprints  of  beams  1  and  2  are  assumed  to  be  1.8  m  apart,  then  an  over¬ 
estimate  of  the  wave  height  by  beam  2  of  7  cm  will  result  in  a  measured  slope  in  the  x 


direction  of 


tan"‘(— )  =  2 

1.8 


O 
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Slope  Time  Series 

Slope  in  y  diracton  Slope  in  xdireclion 


---t  *Oms 


Figure  2.9.3.  Slope  time  series  calculated  with  a  finite  sampling  interval  between 
beam  measurements  vs.  no  sampling  interval  between  beam  measurements.  All  time 
series  were  calculated  from  the  plane  formed  by  the  footprints  along  the  slant  range. 
Incidence  angle,  47*;  Antenna  height,  20  m;  8  =  0*,  Wave  fiequency,  0.2  Hz,  Wave 
direction,  0*,  RLD,  0*. 


An  over-estimate  of  the  height  by  beam  3  will  cause  the  plane  to  have  an  even  greater 
slope  in  the  X  direction  and  also  cause  the  plane  to  have  a  positive  y  slope.  On  the 
back  face  of  the  wave  the  opposite  happens.  The  result  of  this  is  that  the  slope  is 
over-estimated  in  both  the  y  and  x  directions  at  all  times.  Evidence  of  this  can  be 
seen  in  Figure  2.9.3. 

In  Chiqiters  3  and  4,  the  effects  of  these  errors  on  the  VSG's  capability  to 
determine  ocean  parameters  such  as  wave  height  PSD,  the  mean  wave  directional 
spectrum,  and  the  directional  width  spectrum  will  be  analyzed. 
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3.  Determination  of  Mean  Wave  Direction  and  the  Elffect  of  Errors 

3.1  Two-Dimensional  Slope  Distribution  and  Best  Fit  Ellipse 

The  distribution  of  wave  slopes  measured  during  the  SAXON  experiment  is 
nearly  normal  in  both  the  radar  look  direction  and  the  cross  direction  (see  Figure 
3.1.1).  Thus,  the  distribution  of  slopes  is  approximately  bivariate  Gaussian. 


Sy(Deg) 


Sx(Deg) 

Fig.  3.1.1  Histogram  of  x  and  y  slope  components  for  data  run  1749  on  1 1/19/90. 

The  pdf  of  a  bivariate  normal  distribution  is  well  known  and  is  given  by 
[Shanmugan  and  Breipohl]: 

f{x,y)  = - ^  ,  expf^g(x,:y)l  Eq.  3.1.1 


uhere 


we  see  that  the  isopsx>bability  contoxiis  form  a  family  of  ellipses.  The  co^cients  in 
Eq.  3.1.4  are: 

B  =  -pa,Oy 

o  o  ,  Eq.  3.1.5 

C  =  o\  D  =  c(l-p2)o2a2 

The  desired  parameters  of  the  ellipse-the  length  of  the  major  axis  (2a),  the  length  of 
the  minor  axis  (2b),  and  the  ellipse  orientation  (9)->can  be  determined  from  the 
knowledge  of  these  four  coefBcients  [Batschelet,  1981].  The  parameters  of  the 
ellipse  are: 


^  2D  2D 

VA  +  C-R  ~\A-t-C  +  R 


0  =  tan"  ( 


-1.  2B 


A-C-R 


) 


Eq.  3.1.6 


>^iiere  R  is  given  by  v( A  -  C)^  +  4B^  . 

Since  the  moments  of  the  slope  distribution  are  readily  obtainable,  it  is  a  simple 
procedure  to  fit  an  ellipse  to  the  slope  distribution.  The  MATLAB  3.5  program 
ellipsefm  in  Appendix  A  was  written  to  determine  the  parameters  given  in  Eq.  3.1.6, 
and  also  calculates  the  parametric  equations  that  are  used  to  draw  the  ellipse. 


3  J  Mean  Wave  Direction  and  Statistics  of  Slope  Distribution 

The  geometry  of  the  Gaussian-fit  ellipse  yields  infonnation  about  the  ocean 
waves  including  their  primary  direction  of  travel  and  their  degree  of  long>crestedness. 
Figure  3.2.1  is  a  two-dimensional  slope  distribution  and  its  corresponding  Gaussian- 
fit  ellipse  calculated  from  data  recorded  during  the  SAXON-FPN  experiment.  The 
variable  y  in  Eq.  3.1.1  represents  the  radar  look  direction,  w^e  x  represents  the  cross 
direction.  When  the  constant  c  in  Eq.  3.1.3  is  set  equal  to  one,  approximately  40%  of 
the  volume  imder  the  bivariate  normal  distribution  will  lie  above  the  ellipse.  The 
particular  distribution  shown  in  Fig.  3.2.1  is  from  data  taken  between  1649-1739  UTC 
on  Nov.  19, 1990. 

The  orientation  of  the  major  axis  indicates  the  mean  direction  of  wave  travel 
during  the  time  data  were  recorded.  In  the  direction  of  mean  wave  travel,  the  variance 
of  the  wave  slopes  is  greater  than  in  the  Erection  along  the  crest  of  the  waves.  Thus, 
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the  data  tends  to  spread  fiiuther  from  the  origin  in  the  direction  of  wave  travel  than  in 
tlw  cross  direction. 
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Fig.  3.2.1.  Slope  Distribution  and  Best  Fit  Ellipse  for  data  run  1749  on  11/19/90. 
Wave  direction  =  122®  or  302®,  Major  Axis  =  24.0®,  Minor  Axis  »  20.0®.  The  top  of 
the  diagram  is  North. 


If  the  ocean  waves  are  long-crested,  then  the  variance  of  wave  slopes  will  be 
much  greater  in  the  direction  of  wave  travel  than  along  the  crest  of  the  waves.  How 
long  crested  the  waves  actually  are  can  be  determined  by  the  ratio  of  the  major  axis 
length  to  the  minor  axis  length.  The  larger  the  ratio,  the  more  long-crested  the  waves 
are.  From  Fig.  3.2.1,  we  see  that  the  length  of  the  major  axis  is  24.0®  and  the  length 


39 


of  the  minor  axis  is  20.0*’.  Clearly,  the  waves  within  the  radar  footprint  during  run 
1749  were  not  long-crested. 

In  addition  to  giving  the  mean  direction  of  wave  travel,  the  Gaussian-fit  ellipse 
oiables  statistics  of  the  slope  distribution  to  be  easily  obtained.  When  c  in  Eq.  3.1.3 
is  equal  to  one,  the  variance  of  the  slope  in  the  x-direcdon  (East)  is  equal  to  the 
horizontal  distance  from  the  center  of  the  ellipse  to  the  farthest  point  in  the  x-direction 
on  the  ellipse.  In  a  similar  manner,  the  vertical  component  from  the  center  of  the 
ellipse  to  the  farthest  point  in  the  y-direction  on  the  ellipse  is  equal  to  the  variance  of 
the  slope  in  the  y-direction  (North). 

Although  the  orientation  of  the  major  axis  gives  the  mean  direction  of  wave 
travel,  note  that  there  is  a  180**  ambiguity  in  that  direction.  For  example,  in  Fig.  3.2.1 
we  do  not  know  if  the  waves  are  traveling  toward  122**  or  if  they  are  traveling  toward 
302«. 


3  J  Resolving  the  Directional  Ambiguity  using  a  Histogram  of  Slopes 

The  180**  directional  ambiguity  in  Fig.  3.2.1  can  be  resolved  Ity  taking  advantage 
of  the  asymmetry  of  ocean  waves.  Since  the  VSG  recorded  data  at  a  constant  rate  of 
time,  the  VSG  recorded  more  data  from  the  back  side  of  the  waves  than  the  front.  By 
calculating  the  directional  histogram  of  the  two-dimensioiud  slope  time  series,  this 
asymmetry  will  show  up  as  a  shift  in  the  centroid  of  the  histogram.  If  an  ellipse  is  fit 
to  the  histogram  of  the  slope  data,  the  origin  of  the  ellipse  will  be  shifted  in  the  mean 
direction  fiom  which  waves  are  traveling. 

This  procedure  is  illustrated  with  the  data  from  Figure  3.2.1  and  is  shown  in 
Figure  3.3.1.  Due  to  the  asymmetry  of  ocean  waves,  the  center  of  the  histogram  has 
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^lifted  from  the  center.  One  can  see  from  Fig.  3.3.1  that  the  mean  direction  from 
wdiich  die  waves  were  traveling  was  301**. 


1 1/19  1749  Dir.  Histognin  fh  with  LSQ  Method 


Fig.  3.3.1.  Directional  Histogram  ami  Ellipse  fit  with  Least  Square  Method  for  run 
1749  on  11/19/90.  The  deviation  of  the  ellipse  center  from  the  origin  resolves  the 
directional  ambiguity. 

The  ellipse  in  Fig.  3.3.1  indicates  that  the  mean  direction  of  wave  travel  during 
the  time  data  were  recorded  was  30P.  This  differs  fixim  the  302”  suggested  firom  Fig. 
3.2.1  because  the  ellipse  was  fit  with  a  different  method.  In  Fig.  3.2.1,  the  data  are 
nearly  bivariate  Gaussian.  The  moments  of  the  two-dimensional  slope  distribution 
are  used  to  fit  an  ellipse  to  the  data.  Clearly,  the  data  in  Fig.  3.3.1  are  not  bivariate 
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Gaussian.  For  this  reason  tte  method  of  l«ist  squares  was  used  to  fit  an  ellipse  to  the 
data  in  Fig.  3.3.1. 

Given  values  of  x  and  y,  we  want  to  minimize 

M  ^ 

F  =  ?  +  Bxiyi  +  Cyf  +I>Xi  +  Eyi  - 1)  Eq.  3.3.1 

1-1 


Taking  the  derivative  of  F  with  respect  to  each  of  the  five  constants,  X y ,  and  setting 
diem  equal  to  zero  results  in 

yj 

+ Evi  -  iX*? )  “  ® 

1-1 

+  ^/  +  Eyt  -  l)(x/y/ )  =  0 

M 

^iy*  ^y^  +Dxft- Eyi  -  )  =  0  Eq.  3.3.2 

<-i 

In  matrix  finrn  this  can  be  rewritten  as 
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The  MATLAB  3.S  program  Isqfitm  was  written  to  calculate  the  constants  A,  B,  C,  D, 
and  E  given  a  set  of  x  and  y  values. 

After  determining  the  values  of  A,  B,  C,  D,  and  E,  the  parameters  of  die  ellipse 
had  to  be  derived  as  a  function  of  the  constants.  Since  F  was  assigned  to  be  -1,  the 
equatirm  of  the  ellipse  is 

Ax^^Bxy+Cy^-¥Dx-^Ey-\  =  0  Eq.  3.3.6 


Figure  3.3.2.  Geometry  used  to  obtain  ellipse  parameters. 


The  first  step  is  to  write  the  equation  of  the  ellipse  in  the  x'-y*  coordinate  system. 
From  Figure  3.3.2,  we  see  that  the  x'-y*  coordinate  system  is  translated  and  rotated. 
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The  teansfimnation  from  ^  x-y  coordinate  system  to  the  jmmed  coordinate  system  is 
given 


x  =  l  cos(P)  +  x'cos(a)  -  y’  sin(a ) 
y  «  /sin(P)  +  x'sin(a)  +  y  cos(a) 


Substituting  the  values  of  x  and  y  into  Eq.  3.3.6  gives: 

A(I^  cos^  (fi)  +  2i!r'oos(a)oosO)  -  2/y'sin(a)cos(P)  +  x'^  cos^  (a) 

-  2x'^'sin(a)cos(a)  +  y’^  sin^  (o)) 

+B(l^  sin(P)cos(P)  +  fr'sin(a)cos(P)  •»-(y'cos(a)cos(P) + i!x'cos(a)sin(P) 
+  x'^  sin(a)cos(a)  +  x'y'cos^  (a)  -/y'sin(a)sin(P)  -  x'^'sin^  (o) 
sin(a)cos(a)) 

+C(/^  sin^  (P)  +  2fr'sin(a)sin(P)  +  2(y'sin(p)cos(a)  +  x'^  sin^  (a) 

+ 2x'>'sin(a)cos(a) + y’^  cos^  (o))  j 

-«-D(/cos(P)  +  x'cos(o)  -y  sin(a)) 

+£(/sin(p)+x'sin(a)+ycos(o))  -1  =0 


Now,  to  get  the  true  ellipse  equation  we  need  to  eliminate  the  xy  terms.  Doing 
so  gives: 


-2Xx'^'sin(a)cos(a)  +  frx'ycos^  (o)  -  fix^’sin^  (a) 
+2Cic'ysin(a)cos(a)  =  0 


Eq.  3.3.9 


Sinqdifying,  \vr  get 


-1 


Eq.  3.3.10 
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for  the  angle  between  the  x-axis  and  the  x'-axis.  Now  the  equation  of  the  ellipse  in 
the  primed  coordinate  system  is 


( cos^  (a)  +  Bsin(a)cos(a)  +  Csin^  (a))x'^ 

+(i4  sin^  (a)  -  Bsin(a)cos(a)  +  Ccos^  (a))/^ 

+(2i4/  cos(a)cos(P)  +  Bl  sin(a)cosO)  +  B/cos(a)sm(P) 

+  2Cl  sin(a)sinO)  +  Dcos(a)  +  E  sin(a))x' 

+{-2Al  sin(a)cos(P)  +  Bl  cos(a)cosO)  -  Bl  sin(a)sin(P) 

+  2Cl  sin(P)cos(a)  -  Dsin(a)  -t-  E  cos(a))>' 

+Al^  cos^  (P)  +  Bl^  sin(p)cos(p)  +  C/^  sin^  (P) 
+D/cos(P)  +  £/sin(p)  -1  =0 


Since  a,  p,  1,  A,  B,  C,  D,  and  E  are  known,  Eq.  3.3.1 1  can  be  written  as 

Gx'2  +/fy2  =  0 


Dividing  throi^  by  GH  and  reairai^ing 

^  ^  .K 

h'^’gH  ^  GH  GH 


Now,  if  we  complete  the  square  and  do  some  more  rearranging  we  get 


-4KGH+HI^  +GJ^ 
AG^H 


+ 


-4KGH+HI^+GJ^ 

AGH^ 


1 


Eq.  3.3.11 


Eq.  3.3.12 


Eq.  3.3.13 


Eq.  3.3.14 
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Thus,  we  see  that  the  parameters  of  interest  are: 


Center  of  Ellipse 


Eq.  3.3.15 


Length  of  Major  Axis 


l-4KGff+IfI^  +GJ^ 

i  4G^ff 


Eq.3  .16 


Length  of  Mmor  Axis  2  J - - -  Eq.  3.3.17 

V  4Gff^ 


The  MATLAB  3.5  program  dirdist2jn  in  Appendix  A  uses  these  values  to  plot 
the  least-square-fit  ellipse  to  the  histogram  of  slope  data.  If  G^H  is  larger  than  Glfi, 
dien  die  length  of  the  nuyor  axis  will  be  givoi  by  Eq.  3.3.17  and  die  lengdi  of  the 
minor  axis  will  be  given  by  Eq.  3.3.16. 


3,4  Effect  of  Inherent  Measurement  Errors  on  ]>etennination  of  Direction 

In  Chiqiter  2  the  inherent  errors  involved  with  measurements  made  by  the  VSG 
were  analyzed.  Now,  we  would  like  to  know  what  affect  those  errors  have  on 
determining  the  mean  wave  direction  using  the  procedure  described  in  sections  3.1  - 
3.3. 

Figure  3.4.1  shows  a  two-dimensional  slope  distribution  for  a  simulated  ocean 
sur&ce.  In  this  case  and  in  the  following  examples  the  waves  are  coming  fixim  30”, 
and  have  a  fiequeimy  and  amplitude  of  0.2  Hz  and  1.5  m  respectively.  This  is  an 
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ideal  case  in  >^ch  the  slope  is  measured  at  a  single  point.  Ideally,  we  would  like  the 
VSG  to  give  the  same  results,  but  due  to  the  inherent  errors  there  will  be  some 
differences.  The  histogram  is  not  shown  because  with  the  simulated  ocean  surface 
time  is  no  asymmetry  to  take  advantage  of. 
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Figure  3.4.1.  Slope  distribution  of  simulated  ocean  surface.  Slope  was  calculated 
tmalytically  from  a  single  point  Wave  frequency,  0.2  Hz,  Wave  height  =  3  m,  8  =  0®, 
True  wave  direction,  30®;  RLD,  0®. 


Figure  3.4.2  shows  the  slope  distribution  for  the  same  ocean  surfrm  as  in  Fig. 
3.4.1.  The  difference  is  that  in  this  case  the  slope  was  calculated  with  three  vertical 
beams  making  simultaneous  measurements.  Thus,  it  illustrates  the  error  due  to  the 
first  derivative  ^iproximation.  We  see  from  Fig.  3.4.2  that  the  derivative 
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ai^noximation  has  the  effect  of  widening  ellipse  (the  ratio  of  the  nuyor  axis  length  to 
minor  axis  length  is  reduced),  albeit  slightly.  This  occurs  because  the  slope  time 
series  evaluated  from  the  plane  formed  by  the  radar  footprints  is  no  longer  a  true 
sinusoid.  Recall  from  section  2.8  that  the  derivative  approximation  has  the  effect  of 
modulating  the  sinusoidal  frinction  by  a  sine  function. 


Slope  Distributioo  and  Gaussian  Fit  Blipse 


Figure  3.4.2.  Slope  distribution  of  simulated  ocean  surface.  Slope  was  calculated 
from  3  vertical  beams  with  simultaneous  measurements.  Wave  frequency,  0.2  Hz, 
Wave  height  =  3  m,  5  =  0®,  True  wave  direction,  30®;  RLD,  0®. 


Figure  3.4.3  shows  vdiat  h^)pens  when  measurements  are  made  simultaneously 
along  a  slant  range.  Notice  that  the  ratio  of  the  major  axis  length  to  the  minor  axis 


length  is  slightly  less  than  in  Fig.  3.4.2.  This  occurs  because  of  the  increased  error  in 
the  derivative  {q)proximation  when  measuring  along  a  slant  range  (recall  section  2.8). 

In  Fig.  3.4.3,  the  orientation  of  the  major  axis  indicates  that  the  waves  are  coming 
from  30**~the  correct  determination.  Thus,  the  error  due  to  slant-range  measurement 
and  the  derivative  ^proximation  has  no  appreciable  affect  on  the  determination  of  the 
mean  direction  of  wave  travel. 


Sl<^  Distributios  and  Gaussian  Rt  Ellipse 


Figure  3.4.3.  Slope  distribution  of  simulated  ocean  surface.  Slope  was  calculated 
aloi%  a  slant  range  with  simultaneous  measurements.  Incidence  angle,  47";  Antenna 
Height,  20  m;  Wave  frequency,  0.2  Hz;  Wave  height,  3  m;  5  =  0";  RLD,  0";  True 
wave  direction,  30". 

Even  though  slant-range  measurement  causes  the  VSG  to  over-estimate  the  slope 
in  both  directions  at  every  point  in  time  (see  section  2.7),  it  has  no  noticeable  effect 
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on  the  determination  of  the  mean  wave  dilution.  The  reason  for  this  is  that  the  ratio 


remains  ^proximately  the  same  at  all  times.  Thus,  the  mean  direction  of  wave 


travel,  vdiich  is  determined  by  tan 
measurement 


-1 


(?> 

J 


,  is  not  affected  by  slant  range 


Slope  Distributiofi  and  Gaussian  Ht  Ellipse 


Figure  3.4.4.  Slope  distribution  of  simulated  ocean  surface.  Slope  was  calculated 
aloi^  a  slant  range  with  non-simultaneous  measurements.  Incidence  angle,  47**; 
Antenna  Height  20  m;  Wave  fiequency,  0.2  Hz;  Wave  height  3  m;  8  =  0®;  RLD,  0®; 
True  wave  direction,  30®. 

The  most  serious  error,  as  far  as  determining  mean  wave  direction  is  concerned, 
is  due  to  a  finite  switching  time  between  feeds.  Figure  3.4.4  shows  the  results  when 
ail  sources  of  error  are  taken  into  account.  The  slope  is  calculated  fit>m  the  plane 
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foimed  by  the  footprints  whose  locations  are  measured  along  a  slant  range,  non- 
simultaneously.  From  Fig.  3.4.4  we  see  that  the  VSG  would  measure  a  linear,  long- 
crested,  sinusoidal  wave  train  that  is  coming  fiom  30**  as  coming  &om  lb".  In  this 
case,  then,  the  finite  time  between  measurements  causes  the  VSG  to  estimate  the 
mean  wave  direction  14**  to  the  left  of  where  the  waves  are  actually  coming  from. 

Figure  3.4.5  shows  the  results  when  the  VSG  measures  the  same  ocean  surface 
with  33  Hz  switching.  This  shows  that  increasing  the  switching  rate  by  a  frictor  of  3.3 
significantly  improves  the  determination  of  the  mean  wave  direction.  If  greater 
accuracy  is  required.  Figure  3.4.6  shows  the  results  if  the  feeds  are  switched  every  3 
ms.  We  see  that  at  this  rate  the  error  is  virtually  nil. 

In  conclusion,  then,  the  mean  wave  direction  of  ocean  waves  can  be  estimated, 
with  180*’  ambiguity,  by  the  orientation  of  the  major  axis  of  an  ellipse  that  is  fit  to  the 
two-dimensional  slope  distribution.  The  directional  ambiguity  can  be  resolved  by 
taking  advantage  of  the  asymmetry  of  ocean  waves.  When  an  ellipse  is  fit  to  the 
histogram  of  the  slope  distribution,  the  center  of  the  ellipse  is  shifted  in  the  direction 
fiom  wdiich  waves  are  coming. 

The  derivative  approximation  causes  the  width  of  the  ellipse  to  increase,  and 
measuring  along  a  slant  range  causes  the  ellipse  to  widen  even  further.  Non- 
simultaneous  measurement,  however,  has  the  most  significant  effect  on  the 
determination  of  the  mean  wave  direction.  It  causes  the  VSG  to  miscalculate  the 
mean  wave  direction  by  a  significant  amoimt—H**  when  IS”  separates  the  RLD  and 
the  wave  direction. 
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Slope  Distribtttioa  enil  CaimiM  Fit  Blipte 


Figure  3.4.5.  Slope  distribution  of  simulated  ocean  sui&ce.  Slope  was  calculated 
along  a  slant  range  with  10  ms  between  each  measurement.  Incidence  angle,  47**; 
Antenna  Height,  20  m;  Wave  frequency,  0.2  Hz;  Wave  height,  3  m;  8  *  0*;  RLD,  0®; 
True  wave  direction.  30®.  siopeDi«ib«io.«rfo«««R.Eiiip.e 


Figure  3.4.6.  Slope  distribution  of  simulate  ocean  surface.  Slope  was  calculated 
altmg  a  slant  range  with  3  ms  between  each  measurement  Incidence  angle,  47®; 
Antenna  Height  20  m;  Wave  frequency,  0.2  Hz;  Wave  height  3  m;  8  -  0®;  RLD,  0®; 
True  wave  direction,  30®. 


4.  Determination  of  Ocean  Spectra  and  Comparison  with  Pitch-and- 
Roll  Buoy 

4.1  Description  of  Pitch-and-RoU  Buoy 

During  the  SAXON-FPN  experiment  a  WAVEC  pitch-and-roU  buoy  was 
operated  by  the  Federal  Maritime  and  Hydrographic  Agency  of  Hamburg,  Germany 
[Plant  and  Alpers,  1991]  (see  Figure  4.1.1).  Like  the  VSG,  buoys  have  their  own 
sources  of  error.  The  horizontal  motion  of  the  buoy  must  be  accounted  for  [Kuik  et 
al.t  1988].  The  diameter  of  the  buoy  limits  the  length  of  ocean  waves  that  can  be 
accurately  measured.  In  addition,  mooring  beneath  the  buoy  may  affect  its  capability 
to  follow  the  water  surface.  In  order  to  better  understand  the  capability  of  directional 
wave  measurement  systems  to  accurately  determine  principal  wave  parameters, 
Allender  and  others  [1989]  studied  the  performance  of  a  variety  of  such  systems  as 
part  of  the  WADIC  project. 


Figure  4.1.1.  WAVEC  buoy  used  during  the  SAXON-FPN  experiment 


One  such  system  studied  during  the  WADIC  project  was  the  WAVEC  buoy.  The 
WAVEC  buoy  simultaneously  determines  three  quantities— heave,  pitch  angle,  and 
roll  angle.  Records  of  these  time  series  are  used  to  determine  spectral  parameters 
such  as  wave  height  spectrum,  mean  direction  spectrum,  and  directional  spread 
spectrum  (see  section  4.2).  According  to  results  from  the  WADIC  project,  the 
WAVEC  buoy  underestimates  the  wave  height  spectrum  by  10%  -  20%  for  high  sea 
states  and  by  20%  -  50%  for  low  sea  states  for  waves  between  0.10  Hz  and  0.2S  Hz. 
These  numbers  are  comparable  with  other  buoy  systems,  and  may  be  due  to  the  biroy 
going  through  the  wave  A^en  the  wavelength  iq>proaches  the  diameter  of  the  buoy. 
The  mean  wave  direction  at  the  spectral  peak  was  gmerally  within  a  few  degrees  of 
the  best  estimate  data  set  (BEDS),  and  the  wave  spread  estimate  was  not  much 
different  the  BEDS  either.  However,  there  were  a  few  problems  with  the  WAVEC 
buoy  capsizing  in  steep  waves,  and  on  one  occasion  it  actually  left  its  mooring.  Thus, 
the  length  of  the  data  record  used  to  evaluate  the  accuracy  of  the  buoy  was  not  very 
long,  and  the  results  ^uld  be  used  accordingly. 

42  Directional  Distribution  of  Ocean  Waves  and  Longnet-Higgins  Approach 
The  pitch-and-roll  buoy  operated  by  the  Federal  Maritime  and  Hydrographic 
Agency  recorded  the  vertical  acceleration,  the  pitch-angle,  and  the  roll-angle  time 
series.  With  this  information,  the  wave  height  spectrum,  the  mean-wave-direction 
spectrum,  and  the  directional  width  spectrum  can  be  calculated  using  the  method 
developed  by  Longuet-Higgins  et  al.,  [1963].  The  directional  distribution  of  energy 
can  be  ntpanded  in  a  Fourier  series  rei^esentation  as; 
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z>(e.<D)  = 


2it 


CO 


«=l 


(<a )  cos(fi6)  +  6„(a))  sin(n0) 


Eq.  4.2.1 


2n 

a„((o)  =  j  D(6,o>)cos(n0)d9 
0 

2n 

*n(o>)=  jD(0,CD)sin(/i0)dB  Eq.  4.2.2 

0 


The  first  four  Fourier  coefBcients  are  related  to  the  co-  and  cioss>q)ectra  of  die 
time  series  of  wave  heights  and  slopes  by  [Kuik  et  al..  1988]: 


Eq.  4.2.3 

bjio)-  . 

^Cgg  {Cxx  +  Cyy  ) 

Eq.  4.2.4 

^XX-Cyy 

+c 

Eq.  4.2.5 

2C„ 

'■'XX 

Eq.  4.2.6 

i^ieie  C  re|n«sents  co-spectra,  Q  represents  quad-q)ectra. 

z  represents  the  wave 

height,  and  x  and  y  represent  orthogonal  slope  conqwnents.  It  is  also  assumed  that  k. 

die  wave  number,  is  given  by 
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Eq.  4.2.7 


»  ’-a 

With  this  information,  the  spectra  of  interest  are: 

Eq.  4.2.8 
Eq.42.9 

Eq.  4.2.10 

where  Q  (/)  =  . 


wave  height  spectrum 
mean  wave  direction  spectrum 

directional  width  ^)ectrum 


Cx,(/) 

®1  Qzx 

V2-2C,(/) 


43  Comparison  of  VSG  and  Pitch-and>Roll  Bnoy  in  Determining  Ocean  l^ecfra 

Since  die  VSG  has  the  cqiability  to  record  time  series  of  wave  heights  and 
mdxigcmal  slope  components,  the  Loi^itet-Higgins  procedure  can  also  be  used  widi 
the  data  fiom  the  VSG.  The  MATLAB  3.5  program  wiiqiecjn  was  written  to 
determine  the  wave  spectrum,  the  mean  wave  direction  spectrum,  and  the  directional 
width  (about  the  mean)  spectrum  using  the  results  fiom  section  4.2. 

The  CO-  and  cross-spectra  fiom  tlte  VSG  data  were  analyzed  »sing  128 
points  padded  with  zeros  so  that  a  256  point  FFT  was  used.  The  data  were  filtered 
using  5  point  decimation,  and  each  window  of  data  overlapped  die  previous  window 
by  32  points. 

Figure  4.3.1  shows  the  wave  height  spectrum  calculated  by  the  VSG  and  that 
calculated  by  the  pHch-and-roU  buoy.  The  data  were  taken  between  1700  - 1730  UTC 
wi  November  19,  1990.  Figure  4.3.2  shows  the  mean  wave  direction  spectrum 
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calculated  by  the  two  instruments,  and  Figure  4.3.3  ^ws  the  angular  q)read  about 
the  mean  as  a  function  of  frequency.  The  angular  qiread  is  qgproximately  equal  to 
the  rms  angular  deviation  fix>m  the  mean  wave  direction  [Longuet-Higgins,  et  al, 
1963].  Further  examples  of  i^)ectral  parameters  caloilated  by  the  VSG  are  compiled 
in^>pendixC. 


Wave  Heiglit  PSD 


Figure  4.3.1.  One-sided  wave  height  power  ^lectral  density  calculated  from  data 
taken  during  1700-1730  UTC  (m  Nov.  19, 1990.  The  solid  line  was  calculated  from 
VSG  data  and  the  dashed  line  rqnesents  data  obtained  from  the  WAVEC  buoy. 

The  one-sided  wave  height  PSD  and  the  mean  wave  direction  spectrum 
calculated  by  the  two  instruments  compare  frivotably  as  can  be  seen  in  the  figures. 
The  angular  spread,  however,  is  generally  20**  -  40**  greater  for  the  VSG  at  all 
finquendes-most  likely  due  to  noise  and  due  somewhat  to  the  error  induced  from 
slant-range  measuremoit  (see  section  4.4). 


57 


400 


Mm  Wm  DindMa 


Figure  4.3.2.  Mean  wave  direction  spectrum  calculated  from  datf*  taken  during  1700> 
1730  UTC  on  Nov.  19, 1990.  The  solid  line  was  calculated  from  VSG  data  and  the 
dashed  line  rqxesents  data  obtained  from  the  WAVEC  buoy. 


Figure  4  J.3.  Directional  width  qiectrum  calculated  from  data  taken  during  1700- 
1730  UTC  on  Nov.  19, 1990.  The  solid  line  was  calculated  from  VSG  data  and  the 
dashed  line  rqxesents  data  obtained  from  the  WAVEC  buoy. 
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4.4  Effect  of  Inherent  Errors  on  VSG*8  Determination  of  Ocean  Spectra 

The  wave  height  PSD,  the  mean  wave  direction  spectrum,  and  the  directional 
widdi  spectrum  for  an  ideal  measurement  system  are  shown  in  Figures  4.4.1-3 
respectively.  The  time  series  of  wave  heights  and  slopes  were  determined 
analytically  at  a  single  point  The  waves  have  an  amplitude  of  1 .5  m,  a  freqwmcy  of 
.1953125  Hz,  (the  frequency  was  chosen  so  that  10  periods  is  exactly  512  points  with 
10  Hz  sampling)  and  are  coming  from  30**.  The  MATLAB  program  \thspec2.m  usn 
the  Ltmguet-Higgins  method  given  the  time  series  of  wave  heights  and  slopes  to 
calculate  the  spectra  of  interest  A  512  point  FFT  was  used  with  no  padding,  no 
overlq)ping  of  data,  no  decimation,  and  no  windowing. 

The  cross-q)ectral  components,  qi2  and  q23  (equivalmt  to  Qxy  and  Qxx, 
teq)ectively,  in  Eq.  4.2.9)  were  set  to  1  and  0,  respectively  whenever  either 
cmnponenfs  magnitude  was  less  dian  10*^.  This  has  the  effect  of  setting  the  mean 
wave  direction  to  0**  and  the  angular  spread  to  0**  i^ien  data  are  not  present  or  exists 
in  small  quantities.  Figures  4.4.1-3  show  the  results  for  an  ideal  deterministic  noise* 
free  ocean  surfree  >^ien  measured  from  a  single  point 

In  Figures  4.4.4-6,  the  wave  spectra  are  shown  for  the  same  simulated  ocean 
mr&ce.  In  these  figures,  however,  the  measurements  were  made  with  three  beams— 
each  one  measuring  vertically.  Thus,  these  figures  illustrate  the  effect  of  the 
derivative  s^^proximation  on  the  determination  of  such  spectra.  The  only  noticeable 
effect  is  that  the  angular  spread  is  now  2**  at  the  fimdamental  firequency. 

In  Figures  4.4.7*9,  the  spectra  ate  calculated  for  die  same  ocean  surfiice.  In  these 
figures,  diot^  the  measurements  are  made  along  a  slant  range.  The  radar  in  this 
simulatitm  is  mounted  at  20  m  and  the  angle  of  incidence  is  47**.  These  figures,  then, 
illustrate  the  error  due  to  slant-range  measurements.  Notice  the  harmonic  diat  tqqiears 
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Wavdicighl  PSD 


Figure  4.4.1.  One-sided  wave  height  PSD  calculated  from  simulated  ocean  surfa 
data.  Wave  Amplitude,  1.5  m;  Wave  Frequency,  0.1953125  Hz;  Wave  direction,  30 


Wm  OncctiQii 


Figure  4.4.2.  Mean  wave  direction  spectrum  calculated  from  simulated  ocean  surfrict 
data.  Wave  Amplitude,  1.5  m;  Wave  Frequency,  0.1953125  Hz;  Wave  direction,  30®. 
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Figure  4.4.3.  Directional  width  spectrum  calculated  from  simulated  ocean  suifrue 
data.  Wave  Amplitude,  1.5  m;  Wave  Frequoicy,  .1953125  Hz;  Wave  direction,  30*. 
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Figure  4.4.4.  One-sided  wave  height  PSD  calculated  from  simulated  ocean  sur£ace 
data.  Wave  Amplitude,  1.5  m;  Wave  Frequency,  .1953125  Hz;  Wave  direction,  30*. 
Measurements  were  taken  simultaneously  along  three  vertical  beams. 
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Wtve  Direction 


Figure  4.4.S.  Mean  >wave  direction  spectrum  calculated  from  simulated  ocean  surface 
data.  Wave  Amplitude,  1.5  m;  Wave  Frequency,  .1953125  Hz;  Wave  direction,  30**. 
Measurements  were  taken  simultaneously  along  three  vertical  beams. 


Figure  4.4.6.  Directional  width  spectrum  calculated  from  simulated  ocean  surfoce 
data.  Wave  Amplitude,  1.5  m;  Wave  Frequency,  .1953125  Hz;  Wave  direction,  30**. 
Measurements  were  taken  simultaneously  along  three  vertical  beams. 
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Figure  4.4.7.  One-sided  PSD  calculated  from  simulated  ocean  surface  data.  Wave 
Amplitude,  1.5  m;  Wave  Frequency,  .1953125  Hz;  Wave  direction,  30**,  Incidence 
angle,  47°;  Antenna  height,  20  m.  Measuranents  were  made  simultaneously  along  a 
slant  range.  wav*  oncctioR 


Figure  4.4.8.  Mean  wave  direction  spectrum  calculated  from  simulated  ocean  surfru:e 
data.  Wave  Amplitude,  1.5  m;  Wave  Frequency,  .1953125  Hz;  Wave  direction,  30°, 
Incidence  angle,  47°;  Antenna  height,  20  m.  Measurements  were  made 
simultaneously  along  a  slant  range. 
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Figure  4.4.9.  Dii^tional  width  spectrum  calculated  from  simulated  ocean  surface 
dm.  Wave  Amplitude,  1.5  m;  Wave  Frequency,  .1953125  Hz;  Wave  direction,  30* 
IncidOTce  angle,  47*;  Antenna  height,  20  m.  Measurements  were  made 
simultaneously  along  a  slant  range. 
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Figure  4.4.10.  One-sided  PSD  calculated  from  simulated  ocean  surface  Hatn  Wave 
Amplitude,  1.5  m;  Wave  Frequency,  .1953125  Hz;  Wave  direction,  30*,  Incidence 

angle,  47®;  Antenna  height,  20  m.  Measurements  were  made  non-simultaneously 
along  a  slant  range. 
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Figure  4.4.11.  Mean  wave  direction  spectrum  calculated  fit>m  simulated  ocean 
surface  data.  Wave  Amplitude,  1.5  m;  Wave  Frequency,  .1953125  Hz;  Wave 
direction,  30°,  Incidence  angle,  47°;  Antenna  height,  20  m.  Measurements  were  made 
non-simultaneously  along  a  slant  range. 
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Figure  4.4.12.  Directional  width  spectrum  calculated  from  simulated  ocean  surface 
data.  Wave  Amplitude,  1.5  m;  Wave  Frequency,  .1953125  Hz;  Wave  direction,  30°, 
Incidence  angle,  47°;  Antenna  height,  20  m.  Measurements  were  made  non- 
simultaneously  along  a  slant  range. 
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in  the  wave  height  PSD.  Just  as  slant-range  measurement  had  no  appreciable  effect 
on  the  determination  of  mean  wave  direction  with  the  ellipse  fit,  it  also  has  no  effect 
using  the  Longuet-Higgins  s^proach.  The  angular  spread,  however,  is  affected  by 
slant  range  measurements.  Figure  4.4.9  shows  that  the  angular  spread  at  the  second 
harmonic  is  »5.S**.  This  suggests  that  in  a  continuous  spectrum,  there  could  be  a 
significant  increase  in  the  angular  spread  due  to  harmonics. 

Figures  4.4.10-12  show  the  ocean  spectra  taking  into  account  all  inherent  errors. 
Measurements  were  made  non-simultaneously  along  a  slant  range.  The  wave  height 
PSD  is  still  determined  accurately,  but  the  mean  wave  direction  spectrum  has 
dramatically  changed.  Similar  to  the  ellipse  fit,  the  mean  wave  direction  is  now 
approximately  16**  for  most  fioquencies.  Thus,  we  have  further  evidence  of  the 
detrimental  effects  diie  to  non-simultaneous  measurements. 

It  is  not  clearly  understood  >^y  spectral  components  appear  at  all  fioquencies  in 
the  mean  wave  direction  spectrum  and  the  angular  width  spectrum.  Although  512 
data  points  are  used  in  the  FFT,  this  does  not  constitute  a  full  10  periods  of  data. 
Since  there  are  99  ms  between  samples  at  the  same  point,  512  samples  takes  50.688 
sec.  However,  10  periods  of  data  at  .1953125  Hz  takes  51.2  sec.  However,  the  fact 
that  we  do  not  have  a  full  10  periods  of  data  does  not  completely  account  for  the 
additional  spectral  components. 
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5.  Conclusion 

We  have  seen  that,  like  any  system,  the  VSG  has  inherent  errors  associated  with 
its  method  of  measurement.  The  errors  inherent  to  the  VSG  include  phase  shifts  in 
time  soies  of  wave  heights  and  slopes  due  to  measuring  along  a  slant  range, 
modulation  of  the  slope  time  series  due  to  the  first  order  approximation  of  the 
derivative,  and  miscalculation  of  the  mean  direction  of  jvave  travel  due  to  finite 
switching  time  between  the  feeds. 

Measuring  along  a  slant  range  induces  phase  shifts  because  the  measurement  site 
has  some  horizontal  motion.  The  maximum  amount  of  phase  error  in  the  wave  height 
and  slope  time  series  occurs  at  the  crest  and  trough  of  the  wave  while.  The  maximum 
amount  of  error  depends  iqx)n  the  wave  steepness,  the  angle  of  incidence,  and  the 
wave  direction.  It  is  worse  for  steep  waves,  large  angles  of  incidence,  and  for  iqpwave 
or  downwave  look  direction.  Slant>range  measurements  cause  the  slope  time  series  to 
be  overestimated  at  every  point  in  time  in  both  directions. 

The  slant-range  measurements  cause  harmonics  to  appear  in  the  fiequency 
domain.  Even  for  relatively  steep  waves  the  second  harmonic  is  small  («-20  dB),  and 
higher  harmonics  are  even  smaller. 

Approximating  the  slope  at  point  by  a  plane  causes  the  slope  time  series  to  be 
modulated  by  a  sine  function.  The  relative  error  depends  upon  the  wave  number,  k, 
and  the  distance  between  the  radar  footprints.  Slant-range  measurements  cause  the 
distance  between  the  fooqnints  to  be  modulated.  The  error  is  a  maximum  near  the 
trotigh  of  a  wave  when  the  distance  between  the  footprints  is  a  maximum.  For  the 
steepest  waves  the  error  is  about  9%,  but  for  most  of  the  energy-laden  waves  the  error 
is  less  than  2%.  The  derivative  approximation  causes  the  VSG  to  underestimate  tlw 
magnitude  of  the  slope  time  series  at  every  point. 


The  finite  time  between  beam  measurements  seems  to  have  the  most  serious 
effect  on  the  slope  time  series.  Due  to  the  finite  time  between  measurements,  there  is 
phase  error  in  the  second  and  third  measurements.  In  the  wave  height  time  series  this 
does  not  have  much  effect  because  the  three  measurements  are  averaged.  However, 
the  slope  time  series  depends  on  an  accurate  measurement  of  each  position. 
Therefore,  any  error  in  single  measuronent  has  serious  effects. 

The  distribution  of  ocean  waves  measured  during  the  SAXON  experiment  is 
nearly  bivariate  Gaussian.  This  enabled  us  to  fit  an  ellipse  to  the  two-dimensional 
slope  distribution.  The  orientation  of  the  major  axis  indicates  the  mean  direction  of 
wave  travel  with  180”  ambiguity.  In  addition,  statistics  of  the  slope  distribution  are 
easUy  obtained  fi»m  the  ellipse  fit  to  the  distributioiL 

Slant-range  measurements  cause  the  ellipse  width  to  increase  slightly,  but  have 
no  effect  on  the  determination  of  the  mean  wave  direction.  The  derivative 
^jjffoximation  also  causes  the  width  of  the  ellipse  to  increase  slightly.  This  source  of 
error  also  has  no  affect  on  the  determination  of  the  mean  wave  direction.  Non- 
simultaneous  measurements,  however,  cause  the  mean  wave  direction  to  be 
miscalculated.  For  the  particular  example  in  section  3.4,  the  VSG  would  have 
measured  the  waves  as  coming  fix)m  16”  ^en  then  they  were  actually  coming  fiom 
30”. 

The  directional  ambiguity  can  be  resolved  by  taking  advantage  of  the  asyimnetry 
of  ocean  waves.  Since  the  back  side  of  ocean  waves  is  longer  than  the  fit)nt,  and 
since  the  VSG  records  data  at  a  constant  rate  of  time,  more  data  is  recorded  along  the 
back  side  of  the  waves.  Therefore,  a  histogram  of  the  slope  distribution  will  aid  in 
resolving  the  directional  ambiguity.  By  fitting  an  ellipse  to  the  histogram  of  slopes 
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using  the  method  of  least  squares,  the  center  of  the  ellipse  will  indicate  the  mean 
direction  fix>m  which  waves  were  traveling. 

The  Longuet-Higgins  method  for  determining  ocean  spectra  has  been  used  for  a 
number  of  years  with  pitch-and-roll  buoy  data.  Since  the  same  parameters  are 
available  fix>m  the  VSG,  the  Longuet-Higgins  approach  can  be  used  with  the  VSG 
data.  This  enabled  a  comparison  to  be  made  between  the  VSG  and  a  pitch-and-roU 
buoy. 

The  wave  height  spectrum  and  the  mean  wave  direction  spectrum  compare 
&vorably  between  the  two  instruments.  However,  the  angular  spread  about  the  n^an 
is  generally  20**  -  40**  greater  at  all  frequencies  for  the  VSG  than  for  the  pitch-and-roU 
buoy.  There  may  be  several  explanations  for  this.  Higher  frequency  waves  have  a 
shorter  wavelength,  and  thus  are  more  suscq>tible  to  changes  in  the  wind  direction. 
In  addition,  since  the  waves  are  shorter  fewer  measurements  are  along  each  wave 
Mdiich  results  in  more  noise.  As  the  wavelength  approaches  the  size  of  the  buoy,  the 
buoy  may  no  longer  follow  the  surfirce  very  closely. 

The  derivative  ^rproximation  causes  the  angular  spread  to  increase  slightly-from 
roughly  0**  to  2**.  Measuring  along  a  slant  range  causes  a  harmonic  to  appear  in  the 
spectral  parameters.  The  angular  spread  at  the  first  harmonic  is  approximately  S.S**. 
Due  to  non-simultaneous  measurements,  the  mean  wave  direction  was  roughly  16** 
instead  of  the  correct  direction  of  30”.  The  angular  spread  at  the  frmdamental 
fivquency  remained  about  2”. 


69 


5.1  Rccommciidatioiis  for  Further  Study 

It  would  seem  that  the  most  important  item  that  could  be  done  to  improve  the 
VSG's  performance  would  be  to  decrease  the  time  between  measurements  from  each 
feed.  This  should  improve  the  accuracy  of  the  determination  of  the  mean  wave 
direction. 

Although  slant-range  measurements  do  not  seem  to  affect  the  calculation  of  the 
mean  wave  direction,  they  do  cause  the  slope  time  series  to  be  overestimated  at  every 
point  in  both  directions.  This  may  be  important  when  analyzing  the  modulation  of 
die  radar  signal  due  to  tilt  modulation.  The  only  solution  (and  its  only  partial)  is  to 
operate  the  radar  at  snudl  angles  of  incidrace. 

It  seems  that  the  first  order  approximation  of  the  derivative  does  a  good  job  at 
estimating  the  slope  at  a  point  within  the  plane  formed  by  the  radar  foo^irints.  It  has 
little  effect  on  the  wave  height  time  series  or  tl»  slope  time  series. 

All  the  analysis  in  this  pfqier  has  been  done  with  a  deterministic  ocean  surfr»e. 
The  next  step  would  be  to  do  the  same  analysis  with  a  random  ocean  sur&ce  and 
noise  added  to  the  system.  This  would  give  more  realistic  exaiiq>les  of  the  inherent 
errors  involved  with  measuring  wav^  neients  and  vector  slopes  with  the  VSG. 

It  would  be  nice  to  have  an  analytical  description  of  the  error  in  the  slope  time 
series  due  to  non-simultaneous  measurements.  However,  it  is  an  extremely 
complicated  process.  A  finite  time  between  switching  causes  the  measurements  of 
fimn  beam  1  and  beam  2  to  have  some  phase  error.  The  amount  of  phase  error 
d^ends  only  on  the  frequency  of  the  wave  and  the  time  between  sampling  (see  Fig. 
2.9.1).  Once  the  phase  error  has  been  determined,  the  error  in  the  height 
measurement  (vdiich  will  be  a  function  of  time)  can  be  determined.  This  requires  that 
the  height  at  each  footprint  be  known  assuming  that  the  heights  were  measured 
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simultaaeously.  The  differeiKe  in  height  at  each  point  can  be  used  to  detennine  the 
equation  of  the  plane  (w^ch  will  be  a  function  of  time),  and  hence  the  eiror  in  the 
slope  time  series  in  any  direction.  Since  all  three  points  are  required  to  determine  the 
slope  time  soies,  it  is  not  mough  to  determine  the  error  in  the  wave  number  at  each 
point  At  this  point  only  empirical  evideiu:e  is  available  to  illustrate  the  error  due  to 
non-simultaneous  measurements. 
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Appendix  A 

fiinction  [ran,  x,  y,  z]  range(theta,delta,betal2,beta23,gamma,r,phi,f,A) 

%  [ran,  x,  y,  z]  range(theta,delta,betal2,beta23,gamma,r,phi,f,A) 

%  determiiMS  the  range  to  the  ocean  surface  along  the  three  beams  taking 
%  into  account  all  possible  angles.  This  program  solves  the  nonlinear 
%  equation  that  arrises  from  finding  the  point  of  intersection  of  the 
%  radar  beam  with  the  ocean  surface. 

% 

%  theta  -  the  angle  the  radar  was  rotated  about  the  x-axis  (angle  of  inc.) 

%  delta  -  the  angle  the  radar  was  rotated  about  the  z*  axis 
%  betal2  -  the  angle  between  beams  1  and  2 
%  beta23  -  the  angle  between  beams  2  and  3 
%  gamma  -  the  angle  between  beams  1  and  3 
%r  -initial  guess  of  length  of  beam  1,2,  and  3  in  meters 
%  [1, 2, 3] 

%  phi  -  the  angle  between  the  RDL  and  the  direction  fiom  v^ch  tl^  waves 
%  are  coming 

% 

%  Since  there  may  be  multiple  solutions  to  the  nonlinear  equation,  it  is  best 
%  to  guess  the  most  accurate  length  of  r,  while  making  sure  that  the  guess  is 
%  shorter  than  the  resulting  solution. 

% 

%  Chris  Evans,  RSL,  The  University  of  Kansas,  May  26, 1994 
% 


theta  =  theta*pi/180;  %  theta  and  delta  are  measured  positive  fiom  CCW 

delta  -  delta*pi/180; 

betal2»betal2*pi/180; 

beta23  =  beta23*pi/180; 

gamma  »  gamma*pi/180; 

phi  =  phi*pi/180; 


h  =  20;  %  height  of  radar  (meters)  above  mean  sea  surfrao 
%A  -  5.0;  %  amplitude  of  ocean  wave  in  meters 

Rll  cos(delta)*sin(betal2)*cos(gamma-pi/2)  +  sin(delta)*sin(betal2)*sin(gamma- 
pi/2); 

R21 »  sin(delta)*cos(theta)*sin(betal2)*cos(gamma-pi/2)  - 
cos(delta)*cos(theta)*sin(lMtal2)*sin(gamma-pi/2)  +  sin(theta)*cos(betal2); 

R31  -  sin(delta)*sin(theta)*sin(betal2)*cos(gamma-pi/2)  - 
cos(deita)*sin(theta)*sin(betal2)*sin(gamma-pi/2)  -  cos(theta)*cos(betal2); 
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R12-0; 

R22  »  sii^theta); 

R32  =  -co^theta); 

R13  =  -sin(delta)*sm(beta23): 

R23  =  cos(delta)*cos(theta)*sin(beta23)  +  sin(theta)*cos(beta23); 

R33  *  cos(delta)*sin(thcta)*sin(beta23)  -  cos(theta)*cos(bcta23); 

k  =  ((2*pi*f)^2)/9.8;  %spatial  frequency  of  waves  (2*pi/lambda) 

rl(l)  =  0;rl(2)  =  r(l); 
r2(l)  =  0;r2(2)  =  i(2); 
r3(l)  =  0;r3(2)  =  r(3); 

for  t  *  1,  %  t  is  in  tenths  of  seconds;  10/f  gives  1  period 

i  =  2; 

^^e  (abs(rl(i>-rl(i-l))  >  le-13)  |  (abs(r3(i)-r3(i-l))  >  le-13),  %  error  must  be  less 
than  this# 

xl(i)  =  rl(i)*Rl  1;  x2(i)  =  0;  x3(i)  =  r3(i)*R13; 

yl(i)  =  rl(i)*R21;  y2(i)  »  r2(i)*R22;  y3(i)  =  r3(i)*R23; 

zl(i)  =  rl(i)*R31;  z2(i)  =  r2(i)*R32;  z3(i)  =  r3(i)*R33; 

zetal(i)  *  2*pi*f»33*3*t/1000  +  k*sin(phi)*xl(i)  +  k*cos(phi)*yl(i); 

zeta2(i)  =  2*pi*f*33*(3*t+l)/100(H-  k*sin(phi)*x2(i)  +  k*cos(phi)V2(i); 

zeta3(i)  =  2^pi*f*33*(3*t+2yi00(H-  k*sin(phi)*x3(i)  +  k*cos(phi)*y3(i); 

seal(i)  =  A*cos(zetal(i))  -  h; 

sea2(i)  =  A*cos(zeta2(i))  •  h; 

sea3(i)  =  A*cos(zeta3(i))  -  h; 

Fl(i)  =  zl(i)-seal(i);  F2(i)  =  z2(i)  -  sea2(i);  F3(i)  =  z3(i)  -  sea3(i); 

FdCTl(i)  *  Ml  +  A*sin(zetal(i))*^*sin(phi)*Rl  1  +  k*cos(phi)*R21); 

Fder2(i)  -  R32  +  A*sin(zeta2(i))*(k*sin(phi)*R12  +  k*cos(phi)*R22); 

Fder3(i)  =  R33  +  A*sin(zeta3(i))*(k*sin(phi)*R13  +  k*cos(phi)*R23); 
rl(i+l)  =  rl(i)  -  Fl(i)/Fderl(i); 
r2(i+l)  =  r2(i)  -  F2(i)/Fder2(i); 
r3(i+l)  =  r3(i)  -  F3(i)/Fder3(i); 
i-i+1; 
end 

ran(t,:)  =  [rl(i),  r2(i),  r3(i)]; 
x(t,:)  =  [xl(i-l),  x2(i-l),  x3(i-l)]; 
y(t,:)  =  [yl(i-l),  y2(i-l),  y3(i-l)]; 
z(t,:)  =  [zl(i-l).z2(i-l).z3(i-l)]; 
end 


75 


fiincticm  [ran,  x,  y,  z]  =  slrange(theta,delta,betal2.beta23,ganima,r.phi,f,A) 

%  [ran,  x,  y,  z]  =  slrange(theta,delta,betal2,beta23,ganmia,r,i^,f,A)  determines  the 
%  range  to  the  ocean  surface  along  the  three  beams  taking  into  account  all 
%  possible  angles.  This  inrogram  solves  the  nonlinear  equation  that  arrises 
%  from  finding  the  point  of  intersection  of  the  radar  beam  with  the  ocean 
%  surface. 

% 

%  theta  '  the  angle  the  radar  was  rotated  about  the  x-axis  (angle  of  inc.) 

%  delta  -  the  angle  the  radar  was  rotated  about  the  z*  axis 
%  betal2  -  the  angle  between  beams  1  and  2 
%  beta23  -  the  angle  between  beams  2  and  3 
%  gamma  -  the  angle  between  beams  1  and  3 
%r  -initial  guess  of  length  of  beam  1,2,  and  3  in  meters 
%  [1.2,3] 

%  phi  -  the  ai^e  between  the  ROL  and  the  direction  fix>m  ^^ch  the  waves 
%  are  coming 

% 

%  Since  there  may  be  multiple  solutions  to  the  nonlinear  equation,  it  is  best 
%  to  guess  the  most  accurate  length  of  r,  udiile  making  sure  that  the  guess  is 
%  shorter  than  the  resulting  solution. 

% 

%  Chris  Evans,  RSL,  The  University  of  Kansas,  March  1, 1994 
% 

dieta  =  theta*pi/180;  %  theta  and  delta  are  measured  positive  from  CCW 

delta  =  delta*pi/180; 

betal2»betal2*pi/180; 

beta23*beta23*pi/180; 

gamma  =  gamma*pi/180; 

jrfii  =  phi*pi/180; 

h  »  20;  %  height  of  radar  (meters)  above  mean  sea  surface 
%A  »  5.0;  %  amplitude  of  ocean  wave  in  meters 

R1 1  -  cos(delta)*sin(betal2)*cos(gamma-pi/2)  +  sin(delta)*sin(betal2)*sin(gamma- 
pi/2); 

R21  sin(delta)*cos(theta)*sin(betal2)*cos(gamma-piy2)  - 
cos(delta)*cos(theta)*sin(betal2)*sin(gamma-piy2)  +  sin(theta)*cos(betal2); 

R31  =  sin(delta)'''sin(theta)*sin(l^l2)*cos(gamma-pi/2)  - 
cos(delta)*sin(theta)*sin(betal2)*sin(gamma-pi^)  -  cos(theta)*cos(betal2); 

R12*0; 

R22  »  sin(theta); 
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R32  » -cos(theta); 

R13  *  -sin(delta)*sm(beta23); 

R23  *  cos(delta)*cos(theta)*sin(beta23)  +  sin(theta)*cos(beta23); 

R33  »  cos(delta)*sin(theta)*'sin(beta23)  -  cos(theta)*cos(beta23); 

k  « ((2*pi*f)''2)/9.8;  %spatial  frequency  of  waves  (2*pi/lambda) 

rl(l)  =  0;rl(2)  =  r(l); 
r2(l)  =  0;r2(2)  =  r(2); 
r3(l)  =  0;r3(2)  =  r(3); 

for  t  =  1 :10/f,  %  t  is  in  tenths  of  seconds;  40/f  gives  4  periods 

i  =  2; 

while  (abs(rl(i)-rl(i-l))  >  le-13)  |  (abs(r3(i)-r3(i-l))  >  le-13),  %  error  must  be  less 
than  this# 

xl(i)  =  rl(i)*Rll;  x2(i)  =  0;  x3(i)  =  r3(i)*R13; 

yl(i)  =  rl(i)*R21;  y2(i)  =  r2(i)*R22;  y3(i)  =  r3(i)*R23; 
zl(i)*rl(i)*R31;  z2(i)  =  r2(i)*R32;  z3(i)  =  r3(i)*R33; 
zetal(i)  ®  2*pi*f*t/10+  k*sin(phi)*xl(i)  +  k*cos(phi)*y  l(i); 
zeta2(i)  =  2*pi*f*t/l(H  k*sin(phi)*x2(i)  +  k*cos(phi)*y2(i); 
zeta3(i)  *  2*pi*f*t/10+  k*sin(phi)'*x3(i)  +  k*cos(phi)*y3(i); 
seal(i)  *  A*cos(zetal(i))  -  h; 
sea2(i)  =  A*cos(zeta2(i))  -  h; 
sea3(i)  *  A*cos(2eta3(i))  -  h; 

Fl(i)  =  zl(i).seal(i);  F2(i)  =  z2(i) - sea2(i);  F3(i)  =  z3(i) - sea3(i); 

FdCTl(i)  =  R3I  +  A*sin(zetaI(i))*^*sin(phi)*Rl  1  +  k*cosO?hi)*R21); 

Fder2(i)  =  R32  +  A*sin(zeta2(i))*(k*sin(phi)*R12  +  k*cos(phi)*R22); 

Fder3(i)  =  R33  +  A*sin(zeta3(i))*(k*sin(phi)*R13  +  k*cos(phi)*R23); 
rl(i+l)  =  rl(i)  -  Fl(i)/Fderl(i); 
r2(i+l)  =  r2(i)  -  F2(i)/Fder2(i); 
r3(i+l)  =  r3(i)  -  F3(i)/Fder3(i); 
i-i+1; 
end 

ian(t,:)=[rl(i).r2(i),r3(i)]; 
x(t,:)  =  [xl(i-l),x2(i-l),x3(i-l)]; 
y(t,:)  =  [yl(i-l),y2(i-l),y3(i-l)]; 
z(t,:)  =  [zl(i-l),  z2(i-l),  z3(i-l)]; 
end 
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function  [ran,  x,  y,  z]  =  vrange(x,  y,  f.  A,  phi) 

%  [ran,  x,  y,  z]  =  vrange(x,  y,  f,  A,  phi)  determines  the 
%  range  to  the  ocean  surface  assuming  the  three  points  of  intersection  are 
%  moving  vertically  and  assuming  the  measurements  are  made  simultaneously. 
% 

%x  :  time  series  of  X  coordinates  from  range.m  [xl,x2,  x3] 

%y  :  time  series  of  y  coordinates  from  range.m  [yl,y2,  y3] 

%  f  :  freqtiency  of  ocean  waves  used  in  range.m 

%  A  :  amplitude  of  ocean  wave  in  meters  used  in  range.m 

%  phi  :  direction  from  which  ocean  waves  are  coming  in  deg.  (from  range.m) 

% 

%  ran :  time  series  of  ranges  [rl,  r2,  r3] 

%  X  :  time  series  of  x  coordinates  in  earth  system  [xl ,  x2,  x3] 

%  y  :  time  series  of  y  coordinates  in  earth  system  [y  1 ,  y2,  y3] 

%  z  :  time  series  of  z  coordinates  in  earth  system  [zl ,  z2,  z3] 

% 

%  Chris  Evans,  RSL,  The  University  of  Kansas,  March  28, 1994 
% 


phi  =  phi*pi/180; 

h  =  20;  %  height  of  radar  (meters)  above  mean  sea  surface 
k  =  ((2*pi*f)^2)/9.8;  ^^ospatial  fi:^uency  of  waves  (2*pi/lambda) 
[r,c]  =  size(x); 

xl  =  mean(x(:,l));  x2  =  mean(x(:,2));  x3  =  mean(x(:,3)); 
yl  =  mcan(y(:,l));  y2  =  mean0^(:,2));  y3  =  mean^(:,3)); 

t=  l:r; 
t-f; 


zl  =  A*cos(2*pi*f*t/10  +  k*sin(phi)*xl  +  k*cos(phi)*yl)  -  h; 
z2  =  A*cos(2*pi*f  t/10  +  k*sin(phi)*x2  +  k*cos(phi)*y2)  -  h; 
z3  =  A*cos(2*pi*ft/10  +  k*sin(phi)*x3  +  k*cos(p’  )*y3)  -  h; 

rani  =  sqrt(xl''2  +  yl'^2  +  zl  .^^2); 
ran2  =  sqrt(x2^2  +  y2''2  +  z2  .'^2); 
ran3  =  sqrt(x3''2  +  y3''2  +  z3  .^2); 

X  =  [xl*ones(r,l),  x2*ones(r,l),  x3*ones(r,l)]; 
y  =  [yl*ones(r,l),  y2*ones(r,l),  y3*ones(r,l)]; 
z*[zl,z2,  z3]; 
ran  ==  [ntnl,  ran2,  ran3]; 
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function  [Sy,Sx] »  slopeab2(x,y,z) 

%  [Sy,Sx]  =  slopeab2(x,y^) 

% 

%  Slopeab2  determines  the  slope  time  series  from  the  plane  given  the  x,  y, 
%  and  z,  coordinates  of  the  points  of  intersection  of  the  radar  beams 
%  with  the  ocean  surface. 

% 

%  X  -  X  coordinate  of  the  point  of  intersection,  (n  x  3) 

%  y  -  y  coordinate  of  the  point  of  intersection,  (n  x  3) 

%  z  *  z  coordinate  of  the  point  of  intersection,  (n  x  3) 

% 

%  Sy  -  time  series  of  slope  in  y  (RLD)  direction  (degrees) 

%  Sx  -  time  series  of  slope  in  x  (cross)  direction  (degrees) 

% 

%  Chris  Evans,  RSL,  The  University  of  Kansas,  March  28, 1994 
% 

%  B  &  C  are  vectors  in  the  plane,  and  N  is  the  normal  to  the  plane 

B  =  [(x(:,l)  -  x(:,2)),  (y(:,l)  -  y(:,2)),  (z(:,l)  -  z(:,2))]; 

C  =  [(x(:,3)  -  x(:,2)),  (y(:,3)  -  y(:,2)),  (z(:,3)  -  z(:,2))]; 

%  Normal  components  are  found  from  the  cross  product  of  B  and  C. 

Nx  =  B(:,2).'*C(:,3)  -  C(:,2).*B(:,3); 

Ny  =  B(:,3).*C(:,l).C(:,3).*B(:,l); 

Nz  =  B(:,1).*C(:,2)  -  C(:,1).*B(;,2); 

%  The  slope  of  the  plane  is  the  negative  inverse  of  the  slope  of  the  line. 

Sy  =  180/pi*atan(-Ny./Nz); 

Sx=  180/pi*atan(-Nx./Nz); 
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function  [Sy,Sx]  =  slopq)oi(ran,theta,delta,betal2,beta23,gamma^,phi,f) 

% 

%  [Sy,Sx]  -  slopepoi(ran,theta,delta,betal2,beta23,ganuna^phi,f) 

%  calculates  the  slope  at  the  point  for  which  the  slope  is  being 
%  iq)proximated  with  slopeab.m. 

% 

%  ran  :  n  x  3  time  series  of  ranges 
%  theta  :  radar  angle  of  incidence  in  degrees 
%  delta  :  angle  of  CCW  rotation  about  tl»  z*  axis 
%  betal2  :  angle  between  beams  1  and  2 
%  beta23  :  angle  between  beams  2  and  3 

%  gamma :  angle  between  the  planes  formed  by  beams  1  &  2  and  2  &  3  in  degrees 
%  A  :  amplitude  of  ocean  waves  in  meters  (1/2  waveheight) 

%  phi  :  direction  from  which  the  waves  are  traveling 
%  f  :  frequency  of  ocean  wave  in  Hz 
% 

%  Sy  :  Slope  in  y  direction  in  degrees 
%  Sx  :  Slope  in  x  direction  in  degrees 
% 

%  Chris  Evans,  RSL,  The  University  of  Kansas,  March  10, 1994 
% 


theta  -  theta*pi/l  80; 
delta  -  delta*pi/180; 
ganuna  gamma*pi/l  80; 
betal2  =  betal2*pi/180; 
beta23  =  beta23*pi/180; 

[r,c]  *  size(ran); 

%  coordinates  in  the  antenna  system 

xla  =  ran(:,l)*sin(betal2)*cos(pi/2-gamma); 

yla  =  ran(:,l)*sin(betal2)*sin(pi/2-gamma); 

zla  =  -ran(:,l)*cos(betal2); 

x2a  =  zeros(r,l); 

y2a  =  zeros(r,l); 

z2a  -  -ran(:,2); 

x3a  =  zeros(r,l); 

y3a  =  ran(:3)*sin(beta23); 

Z3a  *  -ran(:3)*cos(beta23); 

%  rot  transforms  from  the  antenna  to  the  earth  coordinate  system 
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rot » [cos(delta),  -sin(delta),  0; 

cos(theta)'''sin(delta),  cos(delta)*cos(theta),  -sin(theta); 
sin(theta)*sin(delta),  cos(delta)*sin(dieta),  cos(theta)]; 

fort=  l:r. 

Ant  -  [xla(t),  x2a(t),  x3a(t); 
yla(t),  y2a(t),  y3a(t); 
zla(t),  z2a(t).  z3a(t)]; 


E  =  rot*  Ant; 

xle(t) »  E(l.l);  x2e(t)  =  E(U);  x3e(t)  =  E(l,3); 
yle(t) »  E(2,l);  y2e(t)  =  E(2;2);  y3e(t)  =  E(2,3); 
zle(t)  =  E(3.1);  z2e(t)  =  E(3^);  z3e(t)  =  E(3,3); 


%  coordinates  are  now  in  the  earth  system 

A1  =  (yle(t)-y2e(t))*(z3e(t)-z2e(t))  -  (y3e(t)-y2e(t))*(zle(t)-z2e(t)); 
B1  =  (x3e(t)-x2e(t))*(zle(t)-z2e(t))  -  (xle{t)-x2e(t))*(z3e(t)-z2e(t)); 
Cl  =  (xle(t)-x2e(t))*(y3e(t)-y2e(t))  -  (x3<Kt)-x2e(t))*(yle(t)-y2^t)); 
A2*x3e(t)-x2e(t); 

B2  =  y3e(t)-y2c(t); 

C2  =  z3e(t)-z2e(t); 

A3  =  xle(t)-x2e(t); 

B3=yle(t).y2e(t); 

C3  =  zle(t)-z2e(t); 


R  » [Al,  Bl,  Cl;  A2,  B2,  C2;  A3,  B3,  C3]; 

D1  =  Al  *x2e(t)  +  Bl  *y2e(t)  +  Cl  •z2e(t); 

D2  =  0.5*(x3e(t)^2  +  y3e(tr2  +  z3e(tr2  -  x2e(t)^2  -  y2e(tr2  -  z2e(t)^2); 
D3  =  0.5*(xle(tr2  +  yle(tr2  +  zle(tr2  -  x2e(tr2  -  y2e(tr2  -  z2e(t)^2); 

S  =  pi;D2;D3]; 

Q  =  inv(R)*S; 

Px  =  Q(l); 

Py»Q(2); 

k  =  (2*pi*fr2/9.8; 

Sy(t)- 

atan(A*k*cos(iAi)*cos(2*pi*f*t/l  0+k*sin(phi)*Px+k*cos(phi)*Py))*  1 80/pi; 
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► 


Sx(t)  = 

atan(A*k*sin(phi)*cos(2*pi*f*t/10+k*sin(phi)*Px+k*cos(phi)*Py))*  1 80/pi; 
end 

Sy  »  Sy';  Sx  =  Sx'; 
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function  [a,b,theta,Ro,beta,fi«q,ang,X,Y,xr,yT]  =  ellipsef(x,y,dir,binwidth) 
yoELLIPSEF.EUipsefit 

%  [a,b,theta,Ro,beta,freq,ang^Y,xr,yr]  *  ellipsef(x«y,dir,binwidth) 

% 

%  X  -  x-data  vector 

%  y  -  y-data  vector 

%  dir  -  y  axis  orientation  with  respect  to  North  in  degree  (East  is  90) 

%  binwidth  -  bin  width  for  x-y  phase  histogram  (in  degree) 

% 

%  a  -  length  of  major  axis 

%  b  -  length  of  minor  axis 

%  theta  -  major  axis  orientation  in  radians: 

%  theta(l) -standard  trig,  presentation  (with  respect  to  X  axis) 

%  theta(2)  -  with  respect  to  North 

%  Ro  -  distance  from  ellipse  center  to  the  origin 

%  beta  -  location  angle  of  ellipse  center 

%  freq  -  frequency  counts  of  the  phase  histogram 

%  ang  -  bin  angles  of  the  phase  histogram 

%  X  -  x-coordinates  of  &e  ellipse  for  Earth  m^  orientation 

%  Y  -  y-coordinates  of  the  ellipse  for  Earth  mrq)  orientation 

%  XT  -  rotated  x-data  vector 

%  yr  -  rotated  y-data  vector  (yr  axis  points  North) 

%  Samuel  Haimov,  3/23/93,  RSL,  The  University  of  Kansas 

if  nargin  <  4  binwidth  -  18;  end 
if  nargin  <  3  dir  =  0;  end 

%  Rotate  the  coordinates  relative  to  the  RLD  and  North 
r  *  (36()-dir)*pi/180; 

XT  =  x*cos(r)  -  y*sin(r); 
yr = x*sin(r)  +  y*cos(r); 

%  Find  directional  histogram 

biimum  »  fix(360/binwidth); 
angp  *  angle(xr+j*yr); 

[fi«q,ang] «  hist(angp,binnum); 

%  OUpse  fit  (Bivar.  Normal  Distrb.);  A(xr-mxy'2+2B(xr-mx)(yr-my)+C(yr- 
myy'2»D 

cv»cov([xr(:)yi(:)]); 

inx>miean(xr);  my^eai^);  %  center  of  the  ellipse 
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beta  «  0;  Ro  *  0; 
if  abs(mx)>le-6 1  abs(my)>le-6 

beta=‘atai]2(my4nx);  %  ellipse  center  in  polar 

Ro*sqrt(mx'^2+my^2);  %  coordinates 

end 

vx*cv(l,l);  vy=cv(2,2); 

A“vy ;  B*-cv(2, 1 );  C=vx;  %  ellipse  eq.  polinomial  coeffs. 

D*2*(vx*vy-B'^2);  %  coeff.  in  the  ellipse  polinom. 

R«sqrt((A-Cr244*B'^2); 

a=sqrt(2*D/(A+C-R));  %  major  axis 

b*sqrt(2*D/(A+C+R));  %  minor  axis 

theti^l )=atan(2*B/(A-C-R));  %  ellipse  orientation 

%  also  theta=0.5*atan(2*B/(A-C)); 

%  Calculate  ellipse  parametric  equations 
lrtii«0:pi/128:2*pi; 

X»inx-fa*cos(theta(l  )).*cos(phi)-b*sin(theta(  1  )).*sin(phi); 
Y*my+a*sin(theta(l)).*cos(phi)+b*cos(theta(l)).*sin(phi); 

%  Modify  theta  with  respect  to  North 


if thcta(l)  >= -pi  &theta(l)<  pi/2  theta(2)  =  -theta(l)  + pi/2;  end 
ifdicta(l)>*  pi/2  &  theta(l)  <=  pi  theta(2)  = -theta(l)  +  2.5*pi;  end 


function  [a,b,dietaJlo,beta,freq,ang^,Y^,yr]  =  ellipseh(x,y,dir,binwidth) 

%  ELLIPSEH.  Ellipse  fit  of  the  x-y  data  histogram. 

%  [a,b,thetaJlo,b^fi«q,ang,X,Y^,yr]  -  ellipseh(x,y,dir,binwidth) 

% 

%  X  -  x-data  vector  (e.g.,  slope  in  x-direction) 

%  y  .  y-data  vector  (e.g.,  slope  in  y-diiection) 

%  dir  -  y  axis  orientation  with  respect  to  North  in  degree  (e.g.,  RLD) 

%  binwidth  -  bin  width  for  x-y  phase  histogram  (in  degree) 

% 

%  a  -  length  of  major  axis 

%  b  -  length  of  minor  axis 

%  theta  -  major  axis  orientation  with  respect  to  North  (radians) 

%  theta(l) -standard  trig,  presentation  (with  respect  to  X  axis) 

%  theta(2)  -  with  respect  to  North 

%  Ro  -  distance  fiom  ellipse  center  to  the  origin 

%  beta  -  location  angle  of  ellipse  center 

%  fieq  -  fi«quency  counts  of  the  phase  histogram 

Vo  ang  -  bin  angles  of  the  phase  histogram 

%  X  -  x-coordinates  of  the  ellipse  fm  Earth  map  orientation 

%  Y  -  y-coordinates  of  the  ellii»e  for  Earth  nuq>  orientation 

%  XT  -  rotated  x-component  of  the  histogram 

%  yr  -  rotated  y-con^nent  of  the  histogram  (yr  axis  points  North) 

%  Samuel  Haimov,  3/23/93,  RSL,  The  University  of  Kansas 

%  Rotate  the  coordinates  relative  to  the  RLD  and  North 
r  =  (360-dir)*pi/180; 

XT  *  x*cos(r)  -  y*sin(r); 
yr  *  x*sin(r)  +  y*cos(r); 

%  Find  directional  histogram 

binnum  =  fix(360/binwidth); 
angp  =  angle(xt+j*yr); 

[fieq,ang]  =  hist(angp,binnum); 
xr  «  fieq.*cos(ang); 
yr  *  fi«q.*sin(ang); 

%  Ellipse  fit  (Bivar.  Normal  Distrb.);  A(xr-mx)^2+2B(xr-mx)(yr-my)+C(yr- 
my)^2-D 

cv-cov([xr(:)yr(:)]); 

mx>nnean(xr);  my»miean(yr);  %  center  of  the  ellipse 

beta  *  0;  Ito  *  0; 


if  abs(inx)>lc-6 1  abs(my)>lc-6 

betv4tan2(my4iix);  %  ellipse  center  in  polar 

Ro-sqrt(inx''2+my''2);  %  coordinates 

rad 

vx?>h:v(1,1);  vy-cv(2^); 

A*vy;  B*-cv(2,l);  C*vx;  %  ellipse  eq.  polinomial  coefiGs. 

I>=2*(vx*vy-B''2);  %  fiw  coefif.  in  the  ellipse  polinom, 

R-sqrt((A-Cr2+4*B^2); 
aFsqtt(2*D/(A+C-R));  %  major  axis 

h-sqrt(2*D/(A+C+R));  %  minor  axis 

theta(l)=atan(2*B/(A-C-R));  %  ellipse  orientation 

%  also  theta=0.5*atan(2*B/(A‘C)); 

%  rt^icnlate  ellipse  parametric  equations 
I^*K):pi/128:2*pi; 

X=nnx-»^*cos(theta(l)).*cos(phi>b*sin(theta(l)).*sin(phi); 

Y»my+a*siii(theta(l)).*cos(phi)+b*cos(theta(l)).*sin(phi); 

%  Modify  theta  widi  respect  to  North 

if theta(l)  >=  -pi  &  theta(l)  <  pi/2  thcta(2)  =  -thcta(l)  +  pi/2;  end 
if  theta(l)  >*  pi/2  &  thcta(l)  <=  pi  thcta(2)  *  -theta(l)  +  2*pi+pi/2;  end 


function  [a,b,thetaJlo,beta^Y]-slopdist(date,fiiame^dec,inc^l(l,bwidth,po.go) 
%  SLOPDIST:  Plots  the  data  points  and  the  best  fit  ellipse  for  the  slope  dist. 

% 

%  [a,b,theta,Ro,beta^,Y]  =  slopdist(date3Mune.ndec,incgrld,bwidth,po,go) 

% 

%  date  -  day  of  the  month  for  a  particular  run  i.e.  19 
%  fiuune  -  file  name  of  a  particular  run  i.e.  1047 
%  ndec  -  number  of  points  of  decimation.  Default  is  no  decimation. 

%  inc  -  radar  incidence  angle 
%  rid  -  radar  look  direction 

%  bwidth  -  bin  width  of  phase  histogram  in  degrees  (9  is  default) 

%  po  -  print  option:  'n'  -  nothing  (default)  'm'-create  meta  file 
%  's' create  ps  file  'b' create  meta  &  ps  files  'q'>  create  ps  &  print 

%  'p' create  meta  &  ps  files  &  print 'd' dump  ps  file  to  networic  [ninter 

%  go  -  graphics  option:  'op'  -  portrait,  half  page  (default) 

%  'ol' •  landsciqx,  'of  -  portrait,  full  page 

%  2a  -  length  of  major  axis  in  degrees 

%  2b  -  length  of  minor  axis  in  degrees 

%  theta  >  most  likely  direction  of  wave  propagaton  in  degrees 

%  Ro  -  distance  fiom  origin  to  center  of  ellipse 

%beta  -  location  angle  of  ellipse  center 

*/oX  •  X  coordinates  of  best  fit  ellipse 

%Y  >  y  coordinates  of  best  fit  ellipse 

%  Chris  Evans,  4/13/93,  RSL,  The  University  of  Kansas 

if  nargin  *=  5  bwidth  =  9;  po  *  'n';  end 
if  nar^  =  6  po  =  'n';  end 
if  nargin  =  7  go  =  'op';  end 

mg  -  loadsax(fiiame,'r’,ndec);  %  loads  the  range  data 
s  »  slopeab(mg4nc);  %  calculates  the  x  and  y  slope  components 
s  s  -  oiws(length(s),l)*mean(s); 

s  'S;  %  positive  slope  is  back  &ce  of  the  wave 

%  Ellipse  fit  for  the  slope  distribution 

[a,b,theta4lo,beta,freq,ang,X,Y,xr,yr]  =  eUipsef(s(:,l),s(:,2),rld,bwidth); 
%Gnq[diics 

xr«xr*  180/pi;  yr-yr*  180/pi;  %  Convert  to  degrees 


X»X*180/pi;  Y  =  Y*180/pi; 
a*a*180/pi;  b  =  b*180/pi; 


r  =  min(X):max(X);  %  x-coordinate 

s  =  tan(theta(l  ))*rf  Ro;  %  eqn  of  line  for  major  axis 

u  =  min(Y):max(Y);  %  y-coordinate 

t  =  (u-Roytan(pi/2+theta(  1 ));  %  eqn  of  line  for  minor  axis 

axisCnonnal'); 

axis('square'),  axis([-30  30  -30  30]); 

plot^Y,xr,yr,'.'j,s,'-',t,u,'-');  grid; 

%title(['l  \r  niim2str(date) ' '  num2str(fiiame) '  Slope  Dist.  and  Best  Fit  Ellipse']); 
%ylabelCSy  (Deg)’); 

%xlabelCSx(Deg)'); 

tcxt(.25,.25,['Dir '  sprintfC%3.0f.l80/pi*tIieta(2))  *  deg’],’sc’); 
text(.25,.20,['Major  Axis  = '  sprintfC%2.1f  ^*a)  ’  deg'],’sc’); 
text(.2S,.lS,['Minor  Axis  = '  siMintf('%2.1f  ^•b) '  deg'], 'sc'); 
text(.50..85, •North’, 'sc') 

if  po  =  'n' 
pause 
else 

printps([num2str(fiiame)  'dir'],go,po); 
eml 
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function  [A3.C4),E/]  =  lsqfit(x,y) 

%  LSQFIT:  finds  the  least  square  fit  for  a  2nd  order,  2'D  polynomial. 
%  i.e.  Ax^2  +  Bxy  +  Cy''2  +  Dx  +  Ey  +  F  =  0 
% 

%  function  [A,B,C4>^  =  lsqfit2(x,y) 

% 

%x  •  X  vector  of  real  data  values 
%  y  -  y  vector  of  real  data  values 
% 

%  A3)CJDJEJF  -  coefficients  of  the  best-fit  polynomial 
% 

%  Written  5/12/93  by  Chris  Evans,  RSL,  University  of  Kansas 


x  =  x(:);y  =  y(:); 


V  =  [x,'^2  x.*y  y.^2  x  y]; 

V  =v'*v; 


a  =  inv(V); 
b  =  sum(v)'; 
c  =  a*b; 


A  =  c(l,l);  B  =  c(2,l);  C  =  c(3,l);  D  =  c(4,l);  E  =  c(5,l);  F  =  -1; 
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function  [maj,min,alpha]  -  dirdist2(foame,date,inc,rl(l4idec,bwidth) 

%  DIRDIST2:  plots  the  least  square  fit  for  the  dir.  hist,  of  the  slope  data 
% 

%  function  [major4ninor,alpha]  =  dirdist2(fiiame,date,inc4’l(Midec,bwidth) 
% 

%fiiame  -  filename  of  particular  run  i.e.  1047 
%  date  -  date  of  run  i.e.  19 
%  inc  -  radar  angle  of  incidence 
%  rid  -  radar  look  direction 
%  ndec  -  number  of  decimation  points  (default  is  10) 

%  bwidth  -  width  of  bin  in  deg  for  phase  histogram  (9  is  default) 

% 

%  maj  -  length  of  major  axis 
%  min  -  length  of  minor  axis 

%  alpha  -  most  likely  direction  of  wave  prop,  in  deg  (rel.  to  North) 

% 

%  Written  4/20/93  by  Chris  Evans,  RSL,  University  of  Kansas 
% 


mg  -  loadsax(fiiame,'r',ndec); 
s  =  slopeab(mg,inc); 
s  =  s-ones(length(s),l)*mean(s); 

s  =  -s;  %  positive  slope  is  back  face  of  wave 

%  Rotate  the  coordinates  relative  to  the  RLD  and  North 
r  =  (360-rld)*pi/180; 

XT  =  s(:,l)*co^r)-s(:,2)*sin(r); 
yr  =  s(;,l)*sin(r)+s(:,2)*cos(r); 

%  Find  the  directional  histogram 
binnum  =  fix(360/bwidth); 
angp  =  angle(xrfj*yr); 

[fieq,ang]  =  hist(angp,birmum); 

XT  =  fi«q.*cos(ang); 
yr  =  fi«q.*sin(ang); 


mx  -  mean(xr); 
my  =  mean(yr); 

1 »  sqrt(mx^2  +  my^2);  %  Distance  fiom  origin  to  center  of  ellipse 

beta  -  atan2(my4nx);  %  Angle  between  x-axis  and  ray  to  center  of  ellipse 

[A3,C,D,EJ?]  =  lsqfit(xr,yr); 
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alpha  -  O.S*atan(B/(A-C));  %  angle  between  major  axis  and  pos.  x  axis 

ca  *  cos(alpha);  sa  =  sin(alpha);  cb  =  cos(beta);  sb  =  sin(beta); 

%  Constants  needed  to  calculate  ellipse  parameters 

G  =  A*ca'^2  +  B*sa*ca  +  C*sa'^2; 

H  =  A*sa'^2  -  B*sa*ca  +  C*ca''2; 

I  =  2*A*l*ca*cb+B*l*sin(alpha+beta)+2*C*l*sa*sb+D*ca+E*sa; 

J  =  -2*A*l*sa*clH-B’*l*cos(^pha+beta)*-2*C*l*sa*ca-D*sa+E*ca; 

K  =  A*l'^2*(ca'^2)+B*l'^2*sb*cb+C*l'^2*(sb^2)+D*l*cb+E*l*slH-F; 

h  =  -I/(2*G);  k  =  -J/(2*H);  %  center  coord,  of  ellipse 

numer  =  4*K*G*H  +  H*r2  +  G*r2; 
a  =  sqrt(numer/(4*(G^2)*H));  %  XU  of  major  axis 

b  =  s4tt(numer/(4*G*H''2));  %  XU  of  minor  axis 

theta  =  0:pi/180:2*pi; 

%  Ellipse  parametric  equations 

x  =  l*cb+(h+a*cos(theta))*ca-(k+b*sin(theta))*sa; 
y  =  l*slrKh+a*cos(theta))*sa+(k+b'*sin(theta))*ca; 

%  Determine  the  "true"  major  axis 

if  a>b  maj  =  2*a;  min  =  2*b;  alpha  =  (pi/2-alpha)*  180/pi; 
elseif  aipha>0  alpha  =  (pi-alf^)*  1 80/pi;  maj  =  2*b;  min  =  2*a; 
else  alpha  =  -alpha*  180/pi;  maj  =  2*b;  min  =  2*a; 
end 

%  Normalize  area  to  1 
norm  =  sqrt(pi*maj*min); 

x  =  x/noEm;  y  =  y/norm;  maj  =  maj/noim;  min  =  min/norm; 
xr^xr/norm;  yr  =  yr/norm;  h  =  h/norm;  k  =  k/nonn; 

if  h  <  0  alpha  =  alpha  +  180;  end 

%  Graphics 
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axis('square'),  axis([-l  1  -1  1]); 
plot(xr,yr;*'JUc,'o',x,y);grid; 

%title(n  1/*  num2str((late) ' '  nuin2str(foame) '  Dir.  Histogram  fit  with  LSQ 
Method']); 

%ylabelCRelative  Units'); 

^^abel('Relative  Units'); 
text(.5,.9,'North','sc'); 

text(.25,.25,rDir '  sprintf('%3.0f  .alpha) '  deg'],'sc'); 
text(.2S,.20,[’Major  axis  = '  sprintf('%3.2f.maj)],'sc'); 
text(.2S,.15,['Minor  axis  = '  sprintfr%3.2f.min)],'sc'); 
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function  whspec(name,inc,rld^dec4n^overlap) 

%  WHSPEC(name,inc4idec^ld«m4ioveriap)  calculates  the  waveheight  spectrum 
%  the  directional  spectrum,  and  the  angular  spread  about  the  mean 
%  direction  using  the  Longuet-Higgins  approach.  It  plots  the  results 
%  along  with  the  results  fiom  F.  Ziemer  of  Geekstacht,  Germany  ^dio 
%  used  a  pitch-and-roll  buoy. 

% 

%name  :  4  digit  name  of  file  run.  i.e.  for  191 11047  enter  1047 
%inc  :  radar  incidence  angle 

%rld  :  radar  look  direction  in  deg  (North  is  0,  East  is  90,  etc.) 

%  ndec  :  number  of  decimation  points 
%  m  :  3-vaIued  row  vector  where 
%  m(l)  is  the  length  oftheKth  segment  of  data 

%  m(2)  is  the  number  of  FFT  points  used  (use  a  power  of  2) 

%  m(3)  is  the  sampling  interval  in  sec  (0.1  *ndec) 

%  noverliq} :  m(ll-  point  sections  should  overlap  noverlap  points 
% 

%  For  more  information  about  periodogram  analysis,  see  'Hscrete-Time 
%  Signal  Processing  by  Oppenheim  and  Schafer,  pgs  730-742. 

% 

%  Chris  Evans,  RSL,  The  University  of  Kansas,  June  16, 1994 
% 

clg;  hold  off; 

rid  =  rld*pi/180; 

mg] « loadsax(name,'wt’,ndec); 
vh  =  +  wh(:,2)  +  ^(:,3))y3; 

[s] »  slopeab(mg,inc); 
s  «  s  -  ones(length(s),l)*'mean(s); 

Sx*s(:,l);  Sy  =  s(:,2); 

sx  -  Sy*sin(rld)+Sx*cos(rld);  %  transforms  slopes  to  North  and  East 
sy  =  Sy*cos(rld)-Sx*sin(rld); 

[P]  -  spec(wfa,sx4n,noverlap); 
a»P(:,l);b  =  P(:,2);ql3»P(:,3); 

[Q]  ^  spec(wh,sy4n,noverliq)); 

cl  1 »  Q(:,l);  c  =  Q(:^);  ql2  =  Q(:.3); 

[R]  «  spec(sx,sy,nMioverl{q)); 

c33  *  R(:.l);  c22  =  R(:,2);  c23  =  R(:,3); 
clear  a;  clear  b;  clear  c; 


A  *  ql2  7sqrt(cl  1  .*(c22  +  c33));  %  Equations  are  from  Ziemer 

B  *  ql3  ./sqrt(cl  1  .*(c22  +  c33)); 


mwd  (180/pi)*atan2(ql2,ql3);  %  Mean  wave  direction 

angsp  *  (1 80/pi)*sqrt(2-2*sqrt(A.‘^2  +  B.'^2));  %  Angular  Spread 

mwd  »  -mwd  +  90;  %  Coordinate  Transformation 

a  =  find(mwd<0); 
mwd(a) «  mwd(a)  +  360; 

%  This  loop  keeps  adjacent  points  near 
%  each  otto.  i.e.  350  deg,  10  deg. 

%for  i  *  2:m(l)  +  1 ;  %  would  be  plotted  as  350, 370. 

%  if  (mwd(i)  -  mwd(i-l))  >  235,  mwd(i)  -  mwd(i)  -  360;  end 
%  if (mwd(i)  -  mwd(i-l))  <  -235,  mwd(i)  *  mwd(i)  +  360;  end 
%end 

[n^n] »  size(cl  1); 
ifrem(n,2)  =  0 
f  =  (10/ndec)*(-n:2:n-2)/(2*n); 
select  =  [2:n/2  n/2+2:n]; 
else 

f  =  (10/ndec)*(0:n-iy(2*(n-l)); 
select «  2:n-l; 
end 

df-0.5*abs(f(2)-fl[l)); 
zoom(l)  *  find(f  >  0-df  &  f  <  0+df); 
zoom(2)  *  find(f  >  0.4-df  &  f  <  0.4+df); 
select »  zoQm(l):zoom(2); 

plot(^select),cl  l(select)),grid 
title([num2str(name) '  Waveheight  PSD']); 
xlabdCFrequency'); 
ylabelC[m''2/Hz]'); 

%imnq^[num2str(name)  'a'],'op','b'); 
pause 

plot(^select)4nwd(select)),grid 
title([num2str(name) '  Wave  Direction']); 
xlabelCFiequmicy'); 
ylabelCDegree'); 


9^)rintps([num2str(name)  'a],’op','b’); 
pnise 

plot(f(select),ang^select)),grid 
title([nimi2str(name) '  Angular  Spread']); 
xlabelCFrequency'); 
ylabcKDcgycc'); 

%printps([num2str(name)  'b'],'op','b'); 
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1 1/19  1406  Dir.  Histogrun  fit  with  LSQ  Method 
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1  ]/19  1846  Dir.  Histogram  fit  with  LSQ  Method 
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1 1/20  1839  Dir.  Histogram  fit  with  LSQ  Method 
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